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ABSTRACT 


/  This  thesis  studies  the  estimation  from  interarrival  and 
service  time  data  of  the  mean  number  of  customers  in  service 
at  time  t  for  an  M/G/6^  queue.  Two  situations  are  considered. 
In  one  the  parametric  form  of  the  service  time  distribution 
is  known.  In  the  special  case  in  which  the  service  time 
distribution  is  exponential  the  approximate  bias  and  vari¬ 
ance  of  the  estimate  are  derived  and  simulation  is  used  to 
study  an  approximate  normal  confidence  interval  procedure. 
Simulation  is  also  used  to  illustrate  that  assuming  a  wrong 
parametric  model  can  lead  to  misleading  results.  In  the 
other  situation,  the  parametric  form  of  the  service  time 
distribution  is  unknown  and  the  empirical  distribution  of 
the  service  times  is  used  in  the  estimate  of  the  mean  number 
of  customers  in  service.  In  the  case  in  which  the  customer 
arrival  rate  is  known  the  distribution  of  the  estimate  is 
derived  and  an  approximate  normal  confidence  interval  proce¬ 
dure  is  suggested.  The  use  of  the  bootstrap  and  jackknife 
procedure  to  estimate  variability  and  construct  confidence 
intervals  for  the  estimate  is  also  studied  both  analytically 
and  by  simulation.  _  .  -  ... 
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A.  DESCRIPTION  OF  THE  PROBLEM 


The  application  of  probability  theory  to  a  wide  variety 
of  congestion  problems  has  been  described  in  many  papers  and 
books  [Refs.  1,2,3].  Results  of  queueing  theory  are 
presented  in  terms  of  component  distribution  function  i  and 
stochastic  processes  (renewal,  Poisson,  etc)  that  are  taken 
as  known;  only  rarely  are  issues  addressed  that  arise  when 
actual  data  is  to  be  used  as  a  basis  for  inference  from  the 
models;  however,  see  Cox(  1965)  [Ref.  4]. 

The  concern  of  this  thesis  is  inference  problems  for  a 
particularly  simple  queueing  model,  the  M/G/co  queue.  In 
this  model,  customers  arrive  according  to  a  Poisson  process 
with  rate  )\  and  there  are  an  unlimited  number  of 
independent  servers.  Service  times  for  each  server  are 
independent,  identically  distributed  with  distribution 
function  F.  Let  X(t)  be  the  number  of  customers  being 
served  at  time  t.  It  is  well  known  that  if  there  are  no 
customers  being  served  at  time  0,  then 

P[X(  t  )=n]  =  exp[  -M(  t )  ]  (1.1) 

k.  I 

where 

M(  t)  =  X  £?(  s)ds 

with  F(t)=l-F(t)  [Ref.  2].  Thus  the  distribution  of  X(t)  is 
Poisson  and  is  characterized  by  its  mean  M(t). 

In  this  thesis  we  will  assume  that  the  service  time 
distribution  F(t)  is  unknown  and  must  be  estimated  from 
service  time  data  and  that  the  arrival  process  is  known  to 
be  Poisson,  except  possibly  for  its  rate  .  We  will  study 
the  estimation  of  the  mean  number  of  customers  being  served 
at  time  t ,  M( t ) . 


B.  SCOPE  OF  THE  THESIS 


The  purpose  of  this  thesis  is  to  study  the  estimate  of 
the  mean  number  of  customers  being  served  at  time  t  for  a 
M/G/oO  queue.  This  mean  completely  characterizes  the 
distribution  of  the  number  of  customers  being  served  at  time 
t.  We  will  assume  that  the  service  time  distribution  and 
possibly  the  customer  arrival  rate  are  unknown  and  must  be 
estimated  from  data. 

We  generally  divide  the  estimation  method  into  two  cases 
which  we  shall  call  "parametric  estimation"  and 
"nonparametric  estimation".  In  the  parametric  estimation 
case,  a  particular  probabilistic  model  is  specified  for  the 
service  time  distribution  and  the  parameters  of  the 
distribution  are  estimated.  The  resulting  estimate  of  the 
survivor  function  is  then  used  in  the  estimate  of  the 
expected  number  of  customers  being  served  at  time  t.  In  the 
nonparametric  estimation  method,  the  empirical  survivor 
function  is  used  in  the  estimate  of  the  expected  number  of 
customers. 

In  most  cases,  parametric  assumptions  concerning  the 
service  time  distribution  are  difficult  to  justify.  Hence 
nonparametric  estimation  procedure  may  well  be  preferred  to 
parametric  estimation  when  actual  data  is  used.  However, 
the  nonparair'  ''ric  estimates  can  be  expected  to  be  less 
efficient  than  the  parametric  ones. 

The  thesis  is  organized  as  follows.  In  Chapter  II,  the 
transient  distribution  for  the  number  of  customers  being 
served  at  time  t  for  the  M/G/Oo  model  is  described  and  the 
equilibrium  distribution  as  time  goes  to  infinity  is 
obtained.  In  Chapter  III,  we  study  parametric  estimates  of 
the  mean  number  of  customers  being  served  under  several 
assumptions  for  service  time  distributions.  In  the  special 
case  in  which  the  service  time  distribution  is  exponential 
the  approximate  bias  and  variance  of  the  estimate  are 


[•] 


derived  and  simulation  is  used  to  study  an  approximate 
normal  confidence  interval  procedure.  Parametric  estimates 
for  gamma,  mixed  exponential,  and  lognormal  distributions 
are  also  considered.  Simulation  is  used  to  study  the  effect 
of  assuming  a  wrong  parametric  model.  In  Chapter  IV,  a 
nonparametric  estimate  of  the  mean  number  of  customers  being 
served  is  described.  This  estimate  is  based  on  the 
empirical  distribution  of  the  service  times.  In  the  case  in 
which  the  customer  arrival  rate  is  assumed  known  the 
distribution  of  the  nonparametric  estimate  is  derived  and  an 
asymptotic  normal  confidence  interval  procedure  is 
suggested.  The  jackknife  and  Bootstrap  methods  for 
obtaining  approximate  confidence  intervals  are  also 
described.  The  different  estimators  are  compared  by 
simulation.  Chapter  V  describes  the  simulation  and  gives 
the  results. 

In  summary,  this  thesis  studies  the  use  of  estimates  of 
the  service  time  distribution  to  obtain  estimates  of  the 
mean  number  of  customers  being  served  for  a  M/G/Co  queue. 
Both  parametric  and  nonparametric  estimates  are  considered 
and  compared  by  simulation. 


The  M/G/oo  queueing  model  is  specified  by  the  following 
assumptions.  There  are  infinite  number  of  servers. 
Customers  arrive  for  service  according  to  a  Poisson  process 
with  rate  \  .  Service  times  are  nonnegative  independent 
identically  distributed  random  variables  with  distribution 
function  F.  When  a  customer  arrives,  he  immediately  starts 
service. 

Let  X(t)  represent  the  number  of  customers  in  service  at 


time  t. 


It  is  well  known  that  if  there  are  no  customers 


being  served  at  time  0,  then 


P{X(t)=kj  =  -  —  — 


exp[-«\p(  t)  ] 


(2.  1) 


where  p(t)=  [  1-F(  s )  ]  ds:  that  is,  X(t)  has  a  Poisson 

distribution  with  mean  \p(t)  [Ref.  2].  Taking  the  limit  as 
t  -too  in  equation  2.1,  we  obtain  the  equilibrium 


distribution 


re©  K 

lim  P  (X(  t  )=k  J  =  -  -  exp[  -  1-F(  x)dx]  (2.2) 

t*co  k/ 

Thus,  the  limiting  distribution  of  X(  t)  as  t-too  is  also 

.  oo — 

Poisson  with  mean  m,  where  m=  )o  F(x)dx  is  the  mean  service 
time.  Therefore,  the  distribution  of  the  number  of 
customers  being  served  at  time  t  is  Poisson  with  mean 


M(  t)  =  F(  x)dx 


(2.3) 


Here,  the  distribution  of  the  number  '  of  customers  being 
served  at  time  t  is  characterized  by  value  of  its  mean  M(t). 
The  value  of  M(t)  depends  upon  the  service  time  distribution 
which  is  assumed  unknown  and  must  be  estimated  from  data. 


A.  DESCRIPTION 


In  this  chapter,  it  will  be  assumed  that  the  parametric 
form  of  the  service  time  distribution  is  known.  In  this 
case  the  estimation  of  the  mean  number  of  customers  being 
served  at  time  t,  M(t),  can  be  considered  to  be  a  function 
of  the  parameter  estimates  of  the  distribution.  In 
particular,  the  estimate  of  M(t),  when  a  parametric  form  of 
the  service  time  distribution  is  assumed,  is  denoted  by 
Mp(t),  then 

Mp(t)  =*  C  F(s)ds  (3.  1) 

where  F(t)  is  a  survivor  distribution  of  an  assumed 
parametric  form. 

In  this  chapter  the  rate  of  the  arrival  process  will  be 
assumed  to  be  unknown.  Maximum  likelihood  estimates  of  the 
mean  interarrival  times  and  the  mean  service  times  are  used 
in  the  estimate  of  Mp(t). 

Four  parametric  service  time  distributions  will  be 
considered:  the  exponential,  the  gamma,  the  mixed 
exponential,  and  the  lognormal  distribution.  In  the 
exponential  case,  moment  approximations  are  used  to  assess 
the  bias  of  the  estimate  and  to  develop  a  confidence 
interval  procedure  based  on  asymptotic  normality.  The 
performance  of  the  confidence  interval  procedure  is  assessed 
by  simulation. 

In  the  remaining  three  parametric  models,  simulation  is 
used  to  assess  the  performance  of  the  parametric  estimates. 
Another  source  of  error  in  using  a  parametric  estimator  is 
that  the  wrong  parametric  form  may  be  used.  The  effect  of 
using  the  (wrong)  exponential  model  in  these  cases  is  also 
assessed  by  simulation. 


Each  simulation  has  300  replications;  each  replication 
consists  of  50  independent  service  times  from  the  specified 
distribution,  and  50  independent  interarrival  times  from  an 
exponential  distribution.  The  average  relative  bias  and  the 
average  relative  square  error  of  Mp(t)  are  used  to  evaluate 
the  performance  of  the  parametric  estimation  method.  All 
simulations  were  carried  out  on  an  IBM  3033  computer  at  the 
Naval  Postgraduate  School  using  the  LLRANDOMII,  random 
number  generating  package  [ Ref.  6]  . 

B.  EXPONENTIAL  SERVICE  TIME 

In  this  section  it  will  be  assumed  that  the  service  time 
distribution  is  exponential;  that  is,  F(  t)  =  1-EXP( -t/AA.  ), 
where  11  is  an  unknown  parameter  and  must  be  estimated  from 
the  observed  data.  The  maximum  likelihood  estimate  of  is 

A  i  i  L 

H=-~  51  x,,  where  x-  is  the  service  time  of  the  i  customer. 

We  will  so  assume  that  the  rate  of  the  Poisson  arrival 

process  \  is  unknown  and  must  be  estimated. 

The  interarrival  times  of  the  customers  are  denoted  by  y, 

-  Yi  /  •  •  /  y-n  •  Since  the  arrival  process  is  Poisson  with  rate 

X  ,  the  interarrival  times  are  mutually  independent, 

positive  random  variables  with  the  exponential  distribution 

function  having  mean  4-  .  The  maximum  likelihood  estimate 
a.  *  A 

of  \  is  x =n/  X.  y-  .  For  an  exponential  service  time 

distribution,  an  estimate  of  the  mean  number  of  customers  in 
service  at  time  t  for  an  M/G/to  queue  is 

Mp(  t )  =  \  M  (  l-exp[  -t/w,]  )  (3.2) 

The  estimate  is  a  nonlinear  function  of  the  estimated 

A  A 

parameters,  \  and  -u  .  In  most  cases,  when  estimating  a 
function  of  the  estimated  parameters,  bias  is  created  by  the 
nonlinear  relationship  of  the  estimated  parameters. 

Approximate  formulas  for  the  bias  and  variance  of  Mp(t)  will 
now  be  derived. 


\  A  .  H 

Let  $  be  the  mean  service  time,  (9  =-^  SL  xK  ,  i=l,  2, .  .  ,  n, 

a  ^  | 

and  q(  be  the  mean  interarrival  time,  <*=  —  .2  y.  ,  i=l,2,..,n. 
By  assumption,  X;  and  Y*  are  independent.  The  estimate  of 

A 

Mp(t)  can  be  represented  by  a  function  of  the  parameters  c< 
and  £  as  follows: 

M(S<,p)  =  (l-exp[-t/a)  )  (3.3) 

0< 

There  are  no  simple,  exact  formulas  for  the  mean  and 
variance  of  the  quotient  of  two  random  variables.  However, 
there  are  approximate  formulas  which  are  sometimes  useful. 
The  approximation  can  be  obtained  from  the  partial  Taylor 
series  expansions  of  M(o(,^)  about  the  true  means,  o(  and  (3  . 
The  expansion  is 

M(^,|3)  =  M(o/,p)  +  +  ^-M(V,*)«*-<3) 

.A 

+  i  +  Rn  (3.4) 

Since  we  assumed  that  the  arrival  process  and  the  service 
times  are  independent,  the  covariance  terms  turn  out  to  be 
zero  when  we  take  the  expectation  of  both  sides  of  equation 
3. 4.  Thus,  we  get 

E(M(o<,£)]  =  M(o<',(3)  +  ^  -~M(o<  ,a  )Var(<*)  +  -  -~M(<V ,(&  ) Var(£  ) 
+  (3.5) 


where  R»^  converges  to  zero  at  the  rate 
estimate  is 


The  variance 


Var[M(o<,£)]  =  [ -2-M(Qf/j3  )]  *Var(3<) 


[-j3-M(0(^)J 


:  Var(|3 


of 


+  Rn 


(3.  6) 
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with  tending  to  zero  at  the  rate  .  An  approximate 
bias  term,  denoted  by  (9p(t),  can  be  derived  immediately  from 
the  equation  3.5,  that  is,  £3p(  t  )=E[  M(^  ,  (S  )-E[  M(  ,  £  )  ]  ]  . 
Subtracting  (3p(  t)  from  the  parametric  value  to  correct  the 
bias,  leads  to  the  bias  corrected  estimate  of  Mp(t). 

In  order  to  compare  the  two  estimates,  bias  and 
bias-corrected,  we  define  the  following  notation.  Let  be 
the  fraction  of  bias  of  Mp(t)  against  the  true  value  M(t), 
and  0i  be  the  fraction  of  square  error  of  Mp(t)  against  the 
square  of  the  true  value:  that  is, 

0  =  iito 

Met) 

&.= 

Met) 

where  M(  t)  =  \Cf(  s)ds. 


(3.8) 


A  simulation  experiment  was  performed  to  assess  the 
performance  of  the  estimates.  In  the  i*1'  replication,  50 
exponential  interarrival  times  having  mean  1  and  50 
exponential  service  times  having  mean  2  were  simulated  and 


estimates 


Mp(  t )  =\  M  (  l-exp[ -t/u]  ) 


(3.9) 


Mp(t)  = 


Mp(t) 


-  fMo 


(3. 10) 


where 

Pp(t)  =  i  ±  )Var(/9 ) 

were  computed.  The  estimated  values  of  ^  and  p  were  used 
in  the  variance  formulas.  The  simulation  was  replicated  300 
times  and  the  average  relative  biases 


a,t) 

3oo  i  =  j  Ntk) 


(  3.  11) 


ef,  -  --~zL 

300  a-v  M(t> 

and  the  average  square  error 


(3.  12) 


ea(t>  =  T 

300  M  ft') 


(3.  13) 


e:<t)  =  ?:  [ 

1  3oo 


(3.  14) 


were  computed. 

Results  of  the  simulation  are  presented  in  Figures  3. 1 
and  3.2.  In  Figure  3.1,  the  dotted  line  shows  the  average 


Figure  3.2  Average  Relative  Square  Error  of  Mp(t) 
for  Exponential  with  U  =2  at  t=l.  (  n=50  ,  r=300  ) 

relative  bias,  B((t),  as  a  function  of  t  for  the  original 

A  -  C 

estimate  Mp(t).  The  solid  line  shows  9,(t)  for  the  bias 
'  A  c. 

corrected  estimate  Mp(t).  This  superimposed  figure 

indicates  that  0,C(t)  for  the  bias-corrected  estimate  (with 

solid  line)  is  almost  constant  and  is  small.  The  bias 

estimate  produces  large  negative  value  of  9,(  t)  but  9t(  t) 

approaches  a  limiting  value  as  t  *00  .  Figure  3.2  shows  of 

—  —  c 

the  average  relative  square  error  9*(  t)  and  Q^t)  plotted  as 
a  function  of  time.  The  dotted  line  gives  QJ  t)  and  the 
solid  line  is  0^(t).  It  appears  from  figures  that  the 
estimate  of  bias  described  in  equation  3. 5  does  correct  for 
the  bias.  However,  in  Figure  3.2  the  bias-corrected 
estimate  has  a  slightly  higher  relative  square  error  than 
the  original  estimate.  This  higher  relative  square  error 
could  be  due  to  correlation  between  the  estimate  itself  and 
the  estimate  of  its  bias. 

Simulation  was  used  to  assess  the  performance  of  the 
following  confidence  interval  procedure.  In  the  it'1 


TABLE  I 


COVERAGE  AND  LENGTH  OF  1C0( l-«)%  C. I.  FOR 
THE  ORIGINAL  ESTIMATE  (  N  =  50,  R  =  300  ) 


trial 

68  J 

/ 

O 

80  ? 

/ 

O 

90  J 

/ 

O 

Length 
(  s.d) 

C.  R. 

Length 

(s.d) 

C.  R. 

Length 
(  s.d) 

C.  R. 

1 

0.  2234 

10.  67 
72.  00 

0. 2879 

4.  00 
83.  00 

0.  3695 

1.  00 
90.  00 

( 0. 0415) 

17.  33 

( 0. 0535) 

13.  00 

( 0. 0687) 

5.  00 

2 

0. 2237 

9.  33 
73.  33 

0.  2882 

4.  33 
81.  00 

0. 3699 

1.  33 
86.  33 

( 0. 0396) 

17.  33 

( 0. 0510) 

14.  67 

( 0. 0655) 

12.  33 

3 

0.  2224 

11.  00 
68.  67 

0.  2866 

4.  33 
79.  00 

0.  3678 

1.  33 
87.  33 

( 0. 0410) 

20.  33 

( 0. 0529) 

16.  67 

( 0. 0678) 

11.  33 

4 

0.  2212 

9.  00 
66.  67 

0.  2851 

3.  67 
76.  00 

0.  3659 

1.  00 
83.  33 

(  0.  0402) 

24.  33 

(  0.  0519) 

20.  33 

( 0. 0666) 

15.  67 

5 

0. 2207 

13.  33 
62.  33 

0. 2844 

4.  67 
77.  67 

0.  3651 

0.  67 
88.  67 

(  0.  0409) 

24.  33 

( 0. 0527) 

17.  67 

( 0. 0677) 

10.  67 

6 

0.  2269 

11.  33 
70.  00 

0. 2925 

3.  67 
82.  67 

0.  3754 

0.  33 
89.  33 

( 0. 0433) 

18.  67 

( 0. 0558) 

13.  67 

( 0. 0717) 

10.  33 

7 

0. 2223 

10.  67 
69.  00 

0.  2865 

3.  00 
82.  00 

0.  3677 

0.  33 
89.  67 

(  0.  0375) 

20.  33 

(  0.  0483) 

15.  00 

( 0. 0620) 

10.  00 

8 

0.  2246 

7.  67 
71.  67 

0.  2895 

2.  33 
81.  33 

0.  3715 

0.  33 
87.  67 

( 0. 0425) 

20.  67 

( 0. 0548) 

16.  33 

( 0. 0704) 

12.  00 

9 

0. 2232 

11.  00 
63.  67 

0. 2876 

6.  00 
82.  67 

0.  3692 

0.  67 
83.  67 

{  0.  0423) 

25.  33 

( 0. 0545) 

21.  33 

( 0. 0700) 

15.  67 

10 

0.  2270 

13.  67 
67.  67 

0.  2926 

3.  00 
81.  00 

0.  3755 

0.  33 
87.  33 

( 0. 0411 ) 

18.  67 

( 0. 0530) 

16.  00 

( 0. 0680) 

12.  33 

Average 

0.  2235 

10.  77 
68.  50 

0. 2881 

3.  90 
80.  63 

0.  3697 

0.  74 
87.  33 

( 0. 0410) 

20.  73 

( 0. 0528) 

15.  47 

( 0. 0678) 

11.  93 

replication  of  the  simulation,  50  exponential  interarrival 
times  having  mean  1,  and  50  exponential  service  times  having 

A  /n  r 

mean  2  were  generated,  and  Mp(  t)  and  Mp(t)  were  computed. 
The  approximate  variance  of  Mp(t)  was  computed  for  t=l  using 
equation  3.6  with  Rr,=0.  The  100(1-Q()%  confidence  limits  L 
and  U  were  computed  by 


TABLE  II 


COVERAGE  AND  LENGTH  OF  100(  l-oO%  C.  I.  FOR 
THE  BIAS -CORRECTED  ESTIMATE  (  N=50,  R=300  ) 


68 

/ 

0 

80  ? 

/ 

0 

90  ? 

/ 

0 

L.  L  J.  d  X 

Length 

(s.d) 

C.  R. 

Length 
(  s.  d) 

C.  R. 

Length 
(  s.  d) 

C.  R. 

1 

0.  2234 

14.  67 
69.  00 

0.  2879 

6.  67 
80.  67 

0.  3695 

2.  33 
89.  00 

( 0. 0415) 

16.  33 

( 0. 0535) 

12.  67 

( 0. 0687) 

8.  67 

2 

0.  2237 

13.  67 
69.  67 

0. 2882 

5.  67 
80.  33 

0. 3699 

2.  67 
86.  33 

( 0. 0396) 

16.  67 

( 0. 0510) 

14.  00 

( 0. 0655) 

11.  00 

3 

0.  2224 

16.  00 

66.  00 

0.  2866 

7.  67 
76.  33 

0.  3678 

2.  00 
87.  33 

( 0. 0410) 

18.  00 

( 0. 0529) 

16.  00 

( 0. 0678) 

10.  67 

4 

0.  2212 

14.  33 
62.  33 

0.  2851 

6.  00 
74.  33 

0. 3659 

1.  67 
83.  00 

(  0.  0402) 

23.  33 

( 0. 0519) 

19.  67 

( 0. 0666) 

15.  33 

5 

0.  2207 

17.  67 
59.  67 

0.  2844 

9.  00 
74.  67 

0. 3651 

1.  67 
88.  67 

( 0. 0409) 

22.  67 

( 0. 0527) 

16.  33 

( 0. 0676) 

9.  67 

6 

0.  2269 

18.  67 
64.  00 

0.  2925 

6.  33 
81.  33 

0.  3754 

2.  00 

88.  00 

(  0.  0433) 

17.  33 

( 0. 0558) 

12.  33 

(  0.  0717) 

10.  00 

7 

0. 2223 

14.  67 
67.  00 

0.  2865 

5.  67 
80.  33 

0.  3677 

2.  00 
88.  67 

( 0. 0375) 

18.  33 

(0.  0483) 

14.  00 

( 0. 0620) 

9.  33 

8 

0.  2246 

14.  00 
67.  33 

0. 2895 

4.  67 
80.  00 

0.  3715 

0.  67 
88.  00 

( 0. 0425) 

18.  67 

( 0. 0548) 

15.  33 

( 0. 0704) 

11.  33 

9 

0. 2232 

15.  67 
60.  00 

0. 2876 

9.  00 
70.  67 

0. 3692 

1.  67 
83.  33 

( 0. 0423) 

24.  33 

( 0. 0545) 

20.  33 

( 0. 0670) 

15.  00 

10 

0.  2270 

18.  33 
63.  00 

0. 2926 

9.  33 
75.  33 

0.  3755 

1.  67 
86.  00 

_______ 

( 0. 0411) 

18.  67 

( 0. 0530) 

15.  33 

( 0. 0680) 

12.  33 

Average 

0.  2235 

15.  77 
64.  80 

0.  2881 

7.  00 
77.  40 

0. 3697 

1.  84 
86.  83 

( 0. 0410) 

19.  43 

( 0. 0528) 

15.  60 

( 0. 0678) 

11.  33 

(3.  15) 


(3.  16) 
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where  z  is  the  upper  1-  ^  point  of  the  standard  normal 
distribution.  Tables  I  and  II  show  the  results  of  10 
independent  simulations  for  the  original  and  the 
bias-corrected  estimate.  Each  simulation  was  replicated  30C 
times.  Tables  report  the  average  and  standard  deviation  of 
the  normal  confidence  interval  length;  the  proportion,  p  of 
the  intervals  that  covers  the  true  value;  the  proportion  of 
intervals  that  are  too  high,  ( e. g.  M( t)  <  L  );  and  the 
proportion  of  intervals  that  are  too  low,(e.g.  M(t)  >  U  ). 

Since  the  simulation  replications  are  independent,  it  is 
possible  to  assess  the  uncertainty  of  p.  If  the  confidence 
interval  procedure  is  correct,  then  p  should  be  within 
approximately  ±  2  of  l-o(.  The  coverage  rate  in  the 
tables  indicate  that  the  parametric  estimates  tend  to 

A 

underestimated.  Obviously,  the  distribution  of  Mp(t)  is 
skewed  right.  However,  the  confidence  interval  procedure 
works  well,  regardless  for  both  the  original  and  the 
bias-corrected  estimate.  Both  have  a  variable  coverage 
rate.  The  difference  of  performance  between  two  estimates 
is  not  significant. 

C.  OTHER  SERVICE  TIME  DISTRIBUTIONS 

1.  Mixed  Exponential  Service  Time 

In  this  subsection,  service  times  having  a  mixed 
exponential  parametric  form  will  be  considered.  Customers 
arrive  according  to  a  Poisson  process  with  unknown  rate 
\  which  must  be  estimated.  Customers  are  of  two  types;  with 
probability  P,  ,  a  customer's  service  time  is  exponential 
with  mean  m,  ;  with  probability  P* ,  the  customer's  service 
time  is  exponential  with  mean  M*.  If  T  is  the  service  time 
of  an  arbitrary  customer,  then 


P(T  >  t)  =  P  exp[  -t/u  |  +  Pxexp[ -t/Mx ] 


(3.  17) 


where  P(  +  P a  =  1.  In  this  case,  the  mean  number  of 
customers  being  served  at  time  t  is 


m  p(  t )  =  \  \  p,  a,  Ci-  expc-t/u.'j)  Pi.tu  c  i  -  exp  3. 18) 

We  will  assume  that  a  customer's  type  and  service  time  are 


Figure  3.3  Average  Relative  Bias  of  Mp(t)  for 
Mixed  Exponential  with  u,  =2 ,  w*-=.  75,  P,  =.  2  at  t=l 

(  n=50,  r=300  ) 

A  simulation  experiment  was  done  to  assess  the 
performance  of  M?(t).  In  the  i^  replication,  service  times 
for  50  customers  were  generated,  where  the  proportion  P,  of 
them  have  a  type  1  and  the  proportion  PA  have  a  type  2.  A 
random  number  was  drawn  to  determine  the  type  of  customer. 
If  a  customer  was  of  type  i,  a  service  time  was  generated 
from  the  exponential  distribution  with  mean  u;  .  For  the 
simulation,  the  proportion  ?,  is  assumed  known.  The  50 
independent  interarrival  times  were  also  generated.  In  this 


Figure  3.4  Average  Relative  Square  Error  of  Mp(t)  for 
Mixed  Exponential  with  u, =2,  Ul=.  75,  P,  =.  2  at  t=l 

(  n=50,  r=300  ) 

case  ,  P,  =0.2,  AA,  =2 .  ,  ^^.=0.75,  and  \  =1.  The  estimate  Mp(t) 
was  computed;  the  estimate  of  P*  is  the  proportion  of 
customers  that  were  type  i;  the  estimate  of  the  mean  service 
time  Ux  was  the  mean  service  time  of  type  i  customers;  the 
estimate  of  \  was  the  mean  interarrival  time.  The  estimate 

A 

Me(t)  of  M(t)  assuming  that  the  service  time  distribution  is 
exponential  was  also  computed,  that  is. 


M  (  l-exp( -t/u]  ) 


(3.  19) 


with  u  equal  to  the  mean  service  time  and  \  is  equal  to 
the  mean  interarrival  time.  The  procedure  was  replicated 
300  times.  Results  of  the  simulation  appear  in  Figures  3.  3 
and  3.4.  £,(t)  is  the  average  relative  bias  of  the  estimate 

and  §i(t)  is  the  average  relative  square  error  as  computed 
in  equation  3.  11  and  3.  13.  The  solid  line  is  for  the 

A 

correct  parametric  model  estimate  MP(t).  The  dotted  line 

r  A 

represents  the  exponential  model  estimate  M^(t). 


In  both  figures,  we  compare  the  correct  model 
(represented  by  a  solid  line)  with  the  erroneous  exponential 
model.  As  expected,  the  exponential  model  clearly  does  not 
perform  as  well  as  the  mixed  exponential  model.  In  Figure 
3.3,  the  level  of  bias  for  the  exponential  model  is  very 
high  early  in  time,  but  is  reduced,  and  stabilizes  as  t  ■♦oo  . 
Notice  that  the  estimate  using  the  exponential  model  appears 
too  large.  The  average  relative  square  error  of  the 
estimates  is  shown  in  Figure  3. 4.  Although  the  results  of 
exponential  model  shows  a  slightly  higher  mean  square  error 
than  that  of  the  mixed  exponential  model,  the  results  of 
exponential  model  are  not  too  much  worse.  This  is  not 
surprising,  since  as  t  -*cd  the  estimate  just  depends  on  the 
mean. 

2.  Gamma  Service  Time 

In  this  subsection,  the  parametric  form  of  the 
service  time  distributioi;  is  gamma.  The  probability  density 
function  is 

f(t)  =  -L-  tK'*  e~'3t  (3.20) 

Too 

where  k  and  0  are  strictly  positive  parameters  of  the 
distribution,  and  k  is  further  assumed  to  be  an  integer.  By 
successive  integrations  by  parts,  we  get 

K-|  , 

F(t)  =  1-2.  e  -—3  —  (3.21) 

»  I 

its  mean  and  standard  deviation  are 


Thus,  k  is  the  parameter  that  specifies  the  degree  of 
variability  of  the  service  times  relative  to  the  mean. 

Since  the  arrival  process  is  Poisson  with  rate  \  , 
the  interarrival  times  are  mutually  independent,  positive 


random  variables  with  the  distribution  function 
G( y )=1-EXP( -  \  y ) ,  where  \  must  be  estimated  from 
interarrival  data  y;  ,  i=l  to  n.  Thus,  the  Mp(t)  in  this 
case  is  obtained  by  the  successive  integrations  by  parts  of 
the  survival  function  of  the  gamma  distribution,  that  is 


Figure  3.5  Average  Relative  Bias  of  Mp(t) 
for  Gamma  with  3=1,  k=2  at  t=l.  (n=50,  r=300  ) 

The  performance  of  Mp  (t)  was  evaluated  by  the 
simulation.  In  the  simulation,  the  parameter  k  of  the  gamma 
distribution  was  assumed  to  be  known,  but  the  rate  of  the 
arrival  process  is  unknown  and  is  estimated.  Two  simulation 
cases  were  run.  In  the  i**1  replication  of  the  simulation,  50 
independent  service  times  were  generated,  where  the  first 
simulation  case  used  the  gamma  distribution  having  -3=1  and 
k=2  and  the  second  case  used  the  gamma  distribution  having 
3  =0. 5  and  k=4;  50  independent  interarrival  times  having 
(3=1  were  also  generated.  For  k=2 ,  the  estimate  is 


(3.23) 


Figure  3.7  Average  Relative  Square  Error  of  M  (t) 
for  Gamma  with  0  =1,  k=2  at  t=l.  (  n=50,  r=300  ) 

Figures  3.5,  3.6,  3.7,  and  3.8.  The  dotted  lines  in  the 
figures  show  the  results  of  the  erroneous  exponential  model. 
The  solid  line  represents  the  results  of  the  gamma  model. 
The  figures  clearly  show  that  the  performance  of  the 
exponential  model  is  not  as  good  as  that  of  the  gamma  model 
initially.  However,  both  models  have  the  same  equlibrium 
state  for  t  >  15,  approximately.  Figure  3.5  shows  the 
average  relative  bias  of  the  estimate  in  the  case  of  k=2  and 
Figure  3.  6  shows  the  same  in  the  case  of  k=4.  In  both 
figures,  the  erroneous  exponential  model  has  a  high  level  of 
bias  and  the  bias  of  the  gamma  model  is  almost  constant; 
however,  both  models  have  exactly  the  same  limiting  value  as 
t  "f 00  .  The  average  relative  square  error  appears  in  the 
Figures  3. 7  and  3. 8.  Figure  3.  7  represents  the  average 

A  A 

relative  square  error  of  Mp(t)  and  Me(t).  in  case  of  k=2 , 
and  Figure  3. 8  represents  the  same  in  the  case  of  k=4.  In 
Figure  3.7,  the  exponential  model  has  a  poorer  performance 
than  the  gamma  model,  but  the  difference  is  small.  In 


Figure  3.  8,  the  exponential  model  has  a  large  value  of  mean 
square  error  initially  and  the  level  of  mean  square  error 
associated  with  the  exponential  model  grows  as  the  value  of 
parameter  k  is  increased. 


Figure  3.8  Average  Relative  Square  Error  of  Mp(t) 
for  Gamma  with  0 =.  5,  k=4  at  t=l.  (  n=50,  r=300  ) 

3.  Lognormal  Service  Time 

In  this  subsection,  the  service  time  distribution  is 
assumed  to  be  lognormal.  Let  X  be  a  random  variable,  and 
let  a  new  random  variable  Y  be  defined  as  Y=ln  X.  If  Y  has 
a  normal  distribution,  then  X  is  said  to  have  a  lognormal 
distribution.  The  density  of  a  lognormal  distribution  is 
given  by 


f(X)  =  ______  e 

x  J  iir<f 


*p[  -  in  X  -§  )x 


(  3.  26) 


where  -00<  §  <  co  and  <J  >0.  Set  Y=ln  X  -  §  and  by  the 


integration 

PjX^x]  =  F^(x) 


_  _j _  rM-§ 

'  i^TT  ).oo  eK?  C-  /Vi  <r4]  Ay 

- 


(3. 27) 


where  G^(  •  )  is  the  standard  normal  distribution  function. 
If  the  service  time  has  a  lognormal  distribution,  then  the 
estimated  mean  number  of  customers  being  served  at  time  t. 


lp(t),  is  obtained  by  integration-by-parts 
Mp(t)  =  K  { t[  1-F(  t )  ]  +  j*sF(s)ds} 


(  3.  28) 


Substituting  equation  3.26  for  f(t)  and  equation  3.27  for 
F( t)  in  the  equation  3.  28,  we  obtain 

Mp(t)  =  e*PL§t£akx(~^)}  (3.29) 

where  G^(  )  is  the  standard  normal  distribution  function. 

The  lognormal  distribution  is  positively  skewed  and  the 
level  of  skewness  depends  upon  the  value  of  mean  and 
variance  of  the  distribution.  If  the  value  of  mean  is 
decreased  but  the  variance  is  increased,  then  the  shape  of 
distribution  tends  to  be  more  skewed  and  it  approaches  the 
shape  of  the  exponential  model. 

The  performance  of  the  parametric  estimate  was 
assessed  by  simulation.  In  each  replication  of  the 
simulation,  50  independent  lognormal  service  times  and  50 
independent  exponential  interarrival  times  were  generated. 
The  simulation  generated  two  sets  of  the  service  times.  One 
set  of  service  times  is  from  the  lognormal  distrbution  with 
§=0.193  and  cf^l  ,  and  the  other  set  is  from  the  lognormal 
distribution  with  §  =0.  568  and  <fv= 0.  5.  To  estimate  the  mean 
and  variance  from  the  data,  the  logarithm  of  the  service 


Figure  3.9  Average  Relative  Bias  of  Mp(t)  for 
Lognormal  with  $  =.193,  0^=1  at  t=l.  (  n=50,  r=300  ) 
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Thus  the  estimate  is 
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( 3. 30) 


where  X  is  the  estimate  of  the  arrival  rate.  An  estimate 
based  on  an  erroneous  exponential  model 


Mft(t)  =XM(1  -  exp[ -t/A]  ) 


(3.31) 


K  I  ^ 

was  also  computed,  where  -u  = -x  X.  x-. 

n  \n  ' 


Figure  3.11  Average  Relative  Square  Error  of  Mp(t) 
for  Lognormal  with  §=.193,  6’1  =  1  at  t=l.  (  n=50,  r=300  ) 

The  simulation  was  replicated  300  times.  Figure  3. 8 

shows  the  tendency  of  Q,(t),  the  average  relative  bias  of 
/•» 

Mp(t),  for  the  correct  model  (shown  with  a  solid  line)  and 
the  erroneous  exponential  model  (shown  with  a  dotted  line) 
for  the  simulation  with  §=0.193  and  <J>  =  1  in  Figure  3.9,  and 
with  §  =0.  563  and  <1=0.  5  in  Figure  3.  10. 


Figure  3.12  Average  Relative  Square  Error  of  Mp(t) 
for  Lognormal  with  §=.568,  cf^=.  5  at  t=l.  (  n=50,  r=300  ) 

In  both  figures,  the  average  relative  bias  of  the 
exponential  model  is  also  large  initially.  As  expected,  the 
average  relative  bias  of  the  exponential  model  in  Figure 
3.  10  is  larger  than  the  results  in  Figure  3.  9.  As  t  *co  , 
the  exponential  model  shows  better  performance  and  has  a 
limiting  value.  Figures  3.  11  and  3.  12  show  the  tendency  of 
BJt),  the  average  relative  square  error,  for  the  correct 
model  (shown  with  a  solid)  and  the  erroneous  exponential 
model  (shown  with  a  dotted  line).  For  Figure  3.11  the 
parameters,  §=0.193  and  <fi=l,  are  used  to  generate  data. 
And  for  Figure  3.12  the  parameters,  |  =0.  568  and  <f  =0.5,  are 
used.  The  value  of  average  relative  square  error  of  the 
exponential  model  in  Figure  3.  12  is  also  higher  than  the 
results  in  Figure  3. 11. 

D.  SUMMARY 

The  general  conclusions  of  this  chapter  are  that  the 
parametric  estimation  method  is  a  highly  efficient  for 


obtaining  estimates  of  M(t)  whenever  the  correct  assumption 
for  the  model  is  given.  The  structure  of  the  estimate  of 
M(t)  is  clearly  biased  since  it  is  a  nonlinear  function  of 
the  estimated  parameters  for  the  service  time  distribution 
and  the  customer  arrival  rate.  However,  for  the  service 
time  distributions  considered  the  indications  are  that  the 
bias  is  small.  Hence  the  parametric  estimation  method 
performs  very  well,  whether  or  not  the  estimate  is  corrected 
for  bias,  when  the  correct  parametric  form  is  used.  However 
the  performance  of  the  parametric  estimation  is  very  poor 
when  the  wrong  parametric  model  is  used.  For  instance,  the 
erroneous  exponential  model  often  has  a  high  level  of  bias 
and  mean-squared  error.  Notice  that  the  exponential  model 
converges  to  the  same  limiting  value  as  the  correct  model  as 
t  -»oo  in  all  the  cases  considered.  This  is  because  as  t  ->oc 
all  estimates  use  the  mean  service  time  to  estimate  the 
integral  of  the  servivor  functions. 


'j*  v  v  v.v  pn  ^v; 
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IV.  NONPARAMETRIC  ESTIMATION 

A.  DESCRIPTION 

Nonparametric  methods  are  statistical  techniques  which 
are  applicable  regardless  of  the  form  of  the  distribution 
function  that  the  measurement  comes  from.  In  this  chapter, 
these  techniques  will,  for  the  most  part,  be  based  on  the 
order  statistics. 

Let  ,  ...,x ^  denote  a  random  sample  from  a  CDF  F, 

and  let  su  ,  s^  , .  .  .  ,  s^  denote  that  corresponding  order 
statistics.  Then  the  sample  CDF  is  defined  by 

A  i 

Fn(t)  =  —(number  of  st- less  than  or  equal  to  t) 

~  n  1  l-oo  .  i  3  $(;•> 

For  fixed  time  t,  F„(t)  is  a  statistic  since  it  is  a 

A 

function  of  the  sample.  In  fact,  for  fixed  time  t,  F„(t) 
has  the  same  distribution  as  that  of  the  sample  mean  of  a 
Bernoulli  random  variable.  We  know  by  the  central  limit 

A 

theorem  that  Fn(t)  is  a  asymptotically  normally  distributed 
with  mean  F(t)  and  variance  (  ~  )  F(  t )  [  1-F(  t )  ]  . 

Recall  that  (in  chapter  II)  the  mean  number  of  customers 
being  served  at  time  t,  M(t),  is  a  function  of  the  arrival 
rate  \  and  the  survivor  function  of  the  service  times. 
Hence  a  nonparametric  estimator,  denoted  by  MM(t),  can  be 

\  /v 

represented  by  the  estimated  values  of  \  and  F(t).  The 
estimated  survivor  function  of  service  time  is 

A 

F„(t)  =|  1  if  0  *  t  <  su;> 

----  if  s  ^  t  <  s  for  i=l , 2 , . . , n- 1 

n  co  o-*t> 

if  *  st^ 
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Now,  using  the  fact  that  MN(t)  =  X  F(s)ds,  we  obtain  a 
nonparametric  esimate 


MN(t)  =,  \t 


i  r  '  4-  m-K  .  , 

A  [  +  ----  t] 

“V\  ^  'J>  'h 


Xl-’-Xs  ) 

^  J-\  u  ) 


if  0  $  t  <  s 


if  s  $  t  <  s 


if  t  j  s 


where  X  is  the  estimated  arrival  rate. 


(4.  1) 


Note  that  the 


nonparmetric  estimate  has  a  limiting  value  as  t  -?oo  ,  that 

*  X 

is,  lim  Mf,(t)=\m  where  m  is  the  mean  service  time. 

\-f  oo 

In  this  chapter,  we  will  consider  two  different 


situations. 


In  one  case,  we  will  assume  that  the  arrival 


rate  is  known  but  the  distribution  of  service  times  is 
unknown  and  must  be  estimated.  Based  on  this  assumption, 

A 

M,j(t)  is  expressed  simply  in  terms  of  the  order  statistics 
of  the  service  times  as  follows 


iMd  =  M-^rs  *  t] 


(4.2) 


when  s  <t  <  s,.,  .  In  the  Appendix  A,  we  derive  the 

distribution  of  MK(t)  in  this  case.  Its  mean  and  variance 


E[M„(t)l  =xfF(s)ds 


(4.  3) 


Var[  MkJ(  t ) ]  =  -{  j*sF(ds)  -  [  s F (  d s )  ]  2  +  t2F(t)F(t) 
-  2tF(  t)[  \*sF(ds)}  1  I 


(4.4) 


Thus,  M^t)  is  an  unbiased  estimate  of  M»j(t).  Further,  as 

A 

the  sample  size  n  is  increased,  Mw(t)  is  asymptotically 
normal.  Thus  an  approximate  normal  100(1-<X)%  confidence 

A 

interval  for  Mw(t)  is  given  by 
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Mw(  t)  ±  z  J  Var[M  ( t)  ] 


(4.  5) 


where  z is  the  upper  1-  ^  point  of  standard  normal 
distribution  and  Var[Mg(t)]  is  given  in  Appendix  A.  In  this 
chapter,  we  will  also  study  the  jackknife  and  the  bootstrap 
procedures  for  obtaining  confidence  intervals  for  M^( t)  in 
the  case  in  which  \  is  known.  In  the  second  case,  we  will 
assume  that  the  arrival  rate  X  is  also  unknown  and  must  be 
estimated.  Then,  a  nonparametric  estimate  M*i(t)  is  the 
product  of  two  estimates, 

*  '  i  K  m-K. 

M^t)  =\[-n-Xs  +  — t]  (4.6) 

^  tN  * 

where  \  =n/ and  K  is  the  number  of  service  times  that 
are  less  than  or  equal  to  t.  There  are  no  exact  functional 
forms  for  the  mean  and  variance  of  M^  (t)  in  this  case. 
However,  the  jackknife  and  the  bootstrap  methods  can  be  used 
to  obtain  confidence  intervals.  This  will  be  described 
below. 

B.  JACKKNIFE  ESTIMATION  METHOD 

In  this  section,  we  will  study  the  jackknife  procedure 
for  obtaining  a  confidence  interval  for  M^(t).  The 
jackknife  was  first  introduced  by  Quenouille  (1949)  for  the 
purpose  of  reducing  the  estimate  bias,  and  the  procedure  was 
later  utilized  by  Tukey  (1958),  to  develop  a  general  method 
for  obtaining  approximate  confidence  intervals  [Refs.  7,8]. 

The  basic  idea  of  the  jackknife  estimation  method  i s  to 
assess  the  effect  of  each  of  the  groups  into  which  the  data 
have  been  divided,  not  by  the  results  for  that  group  alone, 
but  rather  through  the  effect  upon  the  body  of  data  that 
results  from  omitting  that  group.  The  two  bases  of  the 
jackknife  are  that  we  make  the  desired  calculation  for  all 
the  data,  and  then,  after  dividing  the  data  into  groups,  we 
make  the  calculations  for  each  of  the  slightly  reduced 


bodies  of  data  obtained  by  leaving  out  just  one  of  the 
groups.  A  special  case  of  Jackknife  estimation  is  called 
the  "complete  jackknife  estimation",  where  the  number  of 

1L 

subgroups  is  n  (the  size  of  sample);  the  i  subgroup  is 
obtained  by  deleting  the  i'tv'  observation;  thus  the  size  of 
each  subgroup  is  n-1  [Ref.  9].  Attention  will  be  restricted 
to  complete  jackknife  estimation  in  this  study. 

A  * 

Let  Mn.((t)  be  the  estimated  mean  number  of  customers 
being  served  at  time  t  on  the  portion  of  the  sample  that 
omits  the  i4^  sample.  Let  M  (t)  be  the  corresponding 
estimator  for  the  entire  sample  and  define  the  i^ 
Pseudo-value  by 


MA(  t)  =  nM(„|(t)  -  (n-l)Mvl(t) 


(4.  7) 


A  A 

The  jackknife  estimate  Mj  (t)  and  an  estimate  S3  of  its 
variance  are  given  by 


a  .  'n  a 

M ,(  t)  -  ~  SI  M;(t) 

J  rr\ 


(4.  8) 


* 

r 


[Mi(t)  -  M^t)]2 


(4.  9) 


Tukey  ( 1958)  proposes  that  the  n  estimated  pseudo  values  be 
treated  as  approximately  independent  and  identically 
distributed  random  variables  [Ref.  9].  Hence,  the  statistic 


77T  ( ryb  -  M  ) 

C  ~  ]l/a 


has  an  approximate  t-distribution  with  n-1 
freedom,  which  leads  to  the  approximate 


(4.  10) 
degrees  of 
100(  !-<*,  )% 


confidence  interval 


where  t  is  the  upper  1-  ~  critical  point  of  the 

t-distribution  with  n-1  degrees  of  freedom.  The  confidence 
interval  given  by  equation  4.  11  is  a  function  of  the 
estimated  variance.  In  the  remainder  of  this  section,  we 
will  describe  several  methods  of  implementing  the  confidence 
interval  procedure.  We  will  also  obtain  an  analytic 
expression  for  the  jackknife  estimate  and  its  variance 
estimate  for  the  case  in  which  the  arrival  rate  \  is  known. 

1.  Jackknife  Estimate  with  Known  Arrival  Rate 

In  this  subsection,  the  arrival  rate  is  asssumed 
known.  In  this  case  the  closed  form  expression  for  the 
jackknife  estimate  and  its  variance  estimate  can  be  derived. 

The  nonparametric  estimate  of  the  mean  number  of 

A 

customers  being  served  at  time  t,  M,j(t),  can  be  expressed  in 
terms  of  the  service  times  as  follows: 


M*(t) 


A 

k 


A 

AA-  K 


X  [  S.S..  +  ----t] 

N  l\  Co  r\ 


(4.  12) 


where  s  are  the  order  statistics  of  independent  and 

identically  distributed  random  quantities  from  the  unknown 

A 

probability  distribution  F,  and  the  variable  K  is  the  number 
of  ' s  which  are  less  than  t.  This  equation  shows 

A 

immediately  that  M^(t)  is  the  linear  function  of  the  order 
statistics  of  service  times. 

The  jackknife  estimate  is  based  on  sequentially 
deleting  point  and  recomputing  the  estimator.  Removing 

point  from  data  set  gives  a  different  empirical 

A  \ 

F„_^  with  mass  -qy  at  Su  > 

corresponding  recomputed  value  of 

In  the  jackknife  process,  the 

value  is 


probability  distribution 

'  •  •  '  '  ^-H) 

the  estimate. 


% 


,  .  .  ,  and  a 


i^  pseudo 


Mi(t)  =f  \t 


\t 

if  i  >  K 

^  SU) 

A 

if  i  $  K 

(4.  13) 
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for  the  fixed  time  t.  Accordingly  the  pseudo  values  ( t) 
have  just  K+l  different  values.  The  jackknife  estimate  is 


\ 

m-  K. 


(  4.  14) 


This  result  is  exactly  the  same  as  the  original  estimate. 
This  is  because  the  estimate  M7(t)  is  unbiased.  In  Appendix 
B,  the  jackknife  variance  estimate  is  derived  as  follows: 

K  A  ^  ^ 

Vart  Mj(  t)  ]  =  2tF,(t)[-;-|y  -  [  J-V]  2 


+  tlFJt)FJt)| 


(4.  15) 


where  Fn(t)  is  the  sample  survivor  function.  Comparing 

(71  A 

equation  4.15  with  equation  4.3  we  see  that  (----)[Var  {M 
(  t ) } ]  =  E[  Var [M 3(  t ) J ]  .  Thus  the  jackknife  variance  estimate 
tends  to  be  conservative  in  the  sense  that  its  expectation 
is  greater  than  the  true  variance  of  M,j(t).  We  will  now 
describe  two  selected  procedures  to  obtain  confidence 
intervals  for  the  jackknife  estimate.  Tukey  suggested  that 
the  statistic  in  equation  4. 10  has  an  approximate 
t-distribution  with  n-1  degrees  of  freedom,  which  leads  to 
the  approximate  two-sided  100(  !-(*)%  confidence  interval 


M7(t)  ±  t,_*  I  Var[Ma(  t)] 


(4.  16) 


for  MNj(t),  where  t  ,.a  is  the  upper  1-  critical  point  of 
the  t-distribution  with  n-1  degrees  of  freedom.  However, 

A 

the  n  estimate  pseudo  values  have  just  K+l  different  values. 
Hence,  another  possible  procedure  is  to  adjust  the  degrees 
of  freedom  of  the  t-distribution,  that  is,  subtract  one  from 

A 

the  number  of  different  pseudo  values  (K+l),  and  use  the 

result  as  the  degrees  of  freedom.  The  length  of  confidence 

\ 

interval  generated  using  the  adjusted  degrees  of  freedom  (K) 
is  slightly  wider  than  that  generated  using  the  usual 


degrees  of  freedom  (n-1)  and  the  coverage  rate  should  be 
increased. 


Jackknife  Estimate  with  Unknown  Arrival  Rate 
In  this  subsection,  it  will  be  assumed  that  the  rate 


of  Poisson  arrival  process  is  unknown  and  must  also  be 
estimated.  The  maximum  likelihood  estimate  is  \=n/^y;, 
where  y-  is  the  interarrival  time  between  i  and  (i-1)  n 
customers.  A  nonparametric  estimate  of  mean  number  of 
customers  being  served  at  time  t  is  given  by 


A  A  I  ^  rr\  —  K 

M„(t)  = 


(4.  17) 


where  S^'s  are  the  order  statistics  of  independent, 
identically  distributed  random  quantities  from  the  unknown 
probability  distribution  F.  It  is  assumed  that  the  S{i;' s  and 

A 

Y*_  s  are  independent.  The  variable  K  is  the  number  of  S  '  s 
which  are  less  than  t.  The  data  consist  of  two  independent 
random  samples, 

S,,Sa . ~  F  and  Y,  ,  Yx . Yn  -  Q 

F  and  Q  being  two  possibly  different  distribution  on  the 
real  line  with  Q,  the  exponential  distribution  with  mean  5-  . 

A 

From  equation  4.17,  the  estimate  M^(t)  is  the  product  of  two 
estimates.  One  is  the  function  of  yx  ,  \  =  and  the 

other  is  the  function  of  s<  ,  H(s)  =  ^  ^  s;+  — ~  t.  there  are 
many  possible  ways  to  perform  a  two-sample  jackknife 
procedure.  We  will  call  one  method  the  "paired  sample 
jackknife"  procedure.  Since  the  size  of  both  samples  is  the 
same,  we  make  the  one  set  of  observations  by  pairing 
respective  observations,  that  is,  (s,,y  ),(  sx  ,  yj  ,...,(  ,  y 

).  As  with  the  one-sample  jackknife,  we  estimate  the  Mftll(t) 

A.  X 

for  all  the  data,  and  then,  we  estimate  Mm_,(t)  based  on  the 
remaining  data  obtained  by  leaving  out  just  the  i*V'  pair. 

xL  A 

Thus  the  i  "  pseudo  value  Mx( t)  is 


Mx(t)  =  nM at.  ( t )  -  (n-l)M^(t) 


(4.  18) 


and  the 
given  by 


Mj(t) 


jackknife  estimate  M3(t) 


* 

2 

V\ 


Mx(t) 


and  variance  estimates  are 


A.  X 

ST  = 


<*>  C  -  '3  A^\ 


J1  A  A 

i.[M.(t)  -  M  (  t )  ]  2 

A^  J 

Based  on  these  statistics,  an  approximate  two-sided 
(  l-°0%  confidence  interval  is  given  by 


100 


M  (t)  ±  t,_*  S, 


(4.  19) 


where  t  is  the  upper  1-  point  of  t-distribution  with 

n-1  degrees  of  freedom.  A  second  method  is  called  the 
"separated  sample  jackknife"  procedure.  Since  we  assumed 
that  the  X;'s  and  ' s  are  independent,  we  can  perform  the 
jackknife  procedure  separately  for  each  sample,  and  then, 
estimates  which  combine  jackknife  estimates  and  the 
jackknife  variance  estimate  can  be  computed. 

A  * 

Let  M3l(t)  be  the  jackknife  estimate  of  \  and  be 

A.  ' 

the  jackknife  variance  estimate  for  \.  Let  M3i(t)  be  the 

f  t  ■“  ^ 

jackknife  estimate  of  F(  s)ds  and  Vs  be  its  jackknife 

variance  estimate.  Then  the  combined  jackknife  estimate  of 
M3c(t)  is 

MJC(  t)  =  MJ((t)-MJa(t)  (4.20) 


and  the  combined  jackknife  variance  estimate  is 


A  * 


SJC=  ^Vs  +  V^M^t)]*  +  Va[MTl(t)] 


(4.21) 


The  approximate  two-sided  100  (  l-o()  %  confidence  interval  is 
given  by 


MT„(t)  ±  t  .  ~  ST,  /J~n 


(4.  22) 


where  t  t_a  is  the  upper  1-^  point  of  t-distribution  with  n-1 
degrees  of  freedom. 

C.  BOOTSTRAP  ESTIMATION  METHOD 

Efron(1979)  introduced  the  bootstrap  method  for 
estimating  the  distribution  of  a  statistic  computed  from 
observations  [Ref.  10].  The  bootstrap  estimate  is  obtained 
by  replacing  the  unknown  distribution  by  the  empirical 
distribution  of  the  data  in  the  definition  of  the 
statistical  function.  In  practice,  the  distribution  of  the 
statistic  is  approximated  by  Monte  Carlo  methods. 

For  convenience,  the  arrival  rate  is  assumed  to  be  known 
and  equal  to  1,  then  the  nonparametric  estimate  M,j(t)  is 
just  a  function  of  service  times.  This  is  a  one-sample 
problem.  The  bootstrap  procedure  is  as  follows: 

1.  Suppose  that  the  data  points  Xi  ,xt  , .  .  .  ,  Xr\  are 
independent  observations  from  the  unknown  distribution 
F.  Then  the  true  estimate  is 

M*(t)  =  S0tF(s)ds  (4.23) 

2.  We  can  estimate  the  distribution  F  by  the  empirical 
probability  distribution  Fn . 

F^  :  mass  on  each  observed  data  point  x- 
i  =  l ,  2  . ,n 

3.  The  bootstrap  estimate  of  M^j(t)  is 

Tb(  t)  =  J*  Fn(s)ds  (4.24) 

* 

To  obtain  an  estimate  of  variability  for  Ma(t),  we  procede 
as  follows 

(1)  Construct  Fn  «  the  empirical  distribution  function, 
as  just  described. 

(2)  Draw  a  bootstrap  sample  x;  ,  x  *  , .  .  .  ,  x  by 

independent  random  sampling  from  Fn . 

Notice  that  we  are  not  getting  a  permutation  distribution 

since  the  values  of  X-  are  selected  with  replacement  from 


the  set  ( x4 ,xx, . . . ,xn  ).  As  a  point  of  comparison,  the 
ordinary  jackknife  can  be  thought  of  as  drawing  samples  of 
size  n-1  without  replacement. 


(3)  Compute  an  estimate  of  Mm  (t 
replication,  Mw(t).  that  is, 
evaluated  for  the  bootstrap  sa: 


of  Mm  (t)  for  each  bootstrap 
that  is,  the  value  of  statistic 
tstrap  sample. 


M'(t)  =  +  vIn'lI(x?-t)lt 


(4.25) 


where  I(x^t)=\l  if  x  (  t 

10  otherwise 


(4)  Do  step  (2)  some  large  number  "B"  times  obtaining 
independent ^  bootstrap  replications  M**'  (t),  M 

Based  on  the  bootstrap  replications,  the  approximate 
estimate  of  M^(t)  and  its  variance  are  obtained  by 


A  i  5 

MB(t)  =  -e 


(4.26) 


Var[  Ma(  t )  ]  =  ---S[M^(t)  -  Mft(t)]* 
B-  \ 


(4.27) 


A  formula  for  the  conditional  variance  of  MR(t)  given  the 

original  sample  data  is  derived  in  Appendix  C.  This 

expression  is  given  by  a  a 

Var[  M p. ( t )  ]  =  { -J-  X x.2  -  [  X  x J  2  +  t2  [  -  ]  -K- 

1  B'  '  ‘  aa  1  rr\  A  1  ^  ^  #  1  AA  '  r\ 


(  4.  28) 


Notice  that  the  expected  value  of  the  conditional  variance 
of  the  bootstrap  estimate  is  approximately  equal  to  the 
variance  of  the  nonparametric  estimate  of  M^jEt)  which  is 
derived  in  Appendix  A. 

So  far  we  have  considered  the  problem,  where  the  arrival 
rate  is  known.  The  bootstarp  methodology  also  applies  if 
the  arrival  rate  is  unknown  and  is  estimated  from 


interarrival  data. 


Suppose  the  data  consist  of  a  random 
sample  X=(X,  ,X^_  ,X^  )  from  unknown  service  time 
distribution  F  and  an  independent  sample  Y=( Y,  ,YX  ,Y^) 
from  the  exponential  interarrival  time  distrbution  G  with 
unknown  parameter  .  One  bootstrap  procedure  to  estimate 
the  expected  number  of  customers  being  served  at  time  t  is 
to  construct  F^  and  Gn  ,  the  empirical  probability 
distribution  corresponding  to  F  and  G.  Bootstrap  samples 
X*  ^  F«  ,  i=l,2,...n,  Y*~  Gn,  j=l,2,...n,  are  independently 
drawn,  an  estimate  of  Mw(t) 


M*(  t)  =  — "l—  [  X  x*  *  n  -  1  I(  x*<  t )  ]  t  }  (  4.  29  ) 

is  calculated.  As  before  there  are  a  large  number  B  of 
bootstrap  replications.  For  this  case,  the  bootstrap 
estimate  of  MM(t)  and  its  variance  are  still  given  by 
equations  4. 26  and  4. 27.  There  appears  to  be  no  closed  form 

A. 

of  the  analytical  variance  of  Ma  ( t)  in  this  case.  Now  we 
will  describe  methods  to  obtain  approximate  confidence 

A 

intervals  for  the  bootstrap  estimate  MB(t). 

1.  The  Percentile  Method 

A  simple  method  for  assigning  approximate  confidence 
intervals  to  the  nonparametric  estimate  M,g(t)  is  as  follows: 
Let 


C(  t) 


(4.  30) 


be  the  cumulative  distribution  function  of  the  bootstrap 
distribution  of  M^(t);  B  is  the  number  of  bootstrap 
replications.  For  a  given  0<o<<0.5,  define 

LU)  =  c“(oO  ,  U(*0  =  c"<  i-°0 

A  A 

Usually  denoted  simply  by  L  and  U.  This  definition  runs 
into  complications  when  we  actually  try  to  compute  quantiles 

A  A 

L  and  U  from  a  set  of  bootstrap  replications.  To  overcome 


these  difficulties,  we  order  the  bootstrap  replications  from 

> 

smallest  to  largest,  obtaining  the  sorted  data  (t),  for 


i=l  to  B. 


Letting 


represent  any  fraction  between  0  and 


1;  take  Q(o<)  to  be  t)  whenever  Q  is  one  of  the 


functions 

be  the  (  B  *X  +  0.  5  f 


,  for  i  =  l  to  B.  Thus  L(o< )  turns  out  to 
M*^(t)  and  U(oO  to  be  the  (B  *  (  1-X) 


+  0.  5  f 


The  percentile  method  consists  of  taking 


[  L(X)  ,  U(oO  ] 


(4.  31) 


as  an  approximate  l-2o(  confidence  interval  for  MQ(t)  since 

A  * 

X=C(L),  l-o<=C(U),  the  percentile  method  interval  consists 

of  the  central  1-2  o<  proportion  of  the  bootstrap 
distribution. 

2.  The  Bias- corrected  Percentile  Method 

Efron(1980)  suggests  the  following  bias  correction 
for  the  percentile  confidence  interval  procedure  [Ref.  11]. 

A 

He  argues  that  if  Ma(t)  is  not  the  median  of  the  bootstrap 
replication  distribution,  then  a  bias  correction  to  the 
percentile  method  is  called  for.  To  be  specific,  define 


z0  =  l"[C(MR(t)) 


where  C(  t )  = 


as  in  equation  4. 30,  and 


(4.  32) 


is  the 


cumulative  distribution  function  for  a  standard  normal 
variate.  The  bias  corrected  percentile  method  consists  of 
taking 


C  '  {  jf,  (  2z0  -  z^)J  ,  c'[£(2z0+  z„) 


(4.33) 


as  an  approximate  1-2 o(  central  confidence  interval  for 

A 

Mo>(t).  Here  z^is  the  upper  point  for  a  standard  normal 
£.<  Zx)  =  l-T(. 


Notice  that  if  M8(t)  is  the  median  of  the  bootstrap 
distribution  then  zo=0  and  equation  4.  33  reduces  to  equation 
4.  31,  the  uncorrected  percentile  interval.  However,  even 
small  differences  of  Pr{M*  (t)  £  MB(t)}  from  0.5  can  make 
equation  4.  33  much  different  from  equation  4. 31. 


The  purpose  of  the  simulation  in  this  chapter  is  to 
assess  the  performance  of  the  nonparametic  estimation 
methods,  the  jackknife  and  the  bootstrap.  Since  the 
estimate  of  M(t),  the  mean  number  of  customers  being  served 
at  time  t,  is  a  function  of  the  customer  arrival  rate  and 
the  integral  of  the  survivor  function  of  the  service  time 
distribution,  two  simulations  cases  are  done.  The  first 
simulation  case  was  performed  to  estimate  Mtsj(t),  the 
nonparametric  estimate  of  M(t),  as  a  function  of  the  service 
times  with  the  arrival  rate  assumed  to  be  known  and  set 
equal  to  1.  For  this  case,  the  jackknife  and  bootstrap 
estimate  of  the  variance  were  derived  in  the  chapter  IV,  and 
compared  with  the  numerical  estimate  obtained  by  the 


simulation. 


The  second  case  assumed  that  the  customer 


arrival  rate  is  also  unknown  and  must  be  estimated  using 
interarrival  times. 

In  each  replication  of  the  simulation  for  case  1,  50 
independent  service  times  from  a  specified  service  time 
distribution  were  generated.  For  the  bootstrap  procedure, 
500  bootstrap  replications  were  performed.  The  simulation 
was  replicated  300  times.  For  the  purposes  of  comparison, 
we  considered  four  types  of  service  time  distributions, 
which  were  the  exponential,  the  mixed  exponential,  the 
gamma,  and  the  lognormal  distribution.  The  arrival  process 
is  known  to  be  Poisson  process  with  known  rate  \  =1.  The 
same  generated  service  times  were  used  for  each  estimation 
procedure  in  a  replication.  This  reduces  the  variability  of 
the  differences  in  performance  between  the  procedures.  All 
programming  was  done  on  IBM  3033  computer  at  the  Naval 
Postgraduate  school  using  the  LLRANDOMII,  random  number 
generating  package  [ Ref.  6]  . 


’ABLE  III 


T 

STATISTICAL  DATA  OF  ESTIMATE  OF  M(T) 
N=50  R=300  ( B=500 ) 


True 

value 

Parametric 

Nonparametric 

Correct 

Errors 

Jack 

Boot 

Exponential 
(  «  =  2) 

0. 7869 

0. 7830 
( 0. 0268) 

- 

0. 7820 
( 0. 0451) 

0. 7819 
(  0.  0452) 

0.  5871 
{ 0. 0023 ) 

0.  6246 
( 0. 0026) 

0.  5991 
( 0. 0512 ) 

0. 5989 
( 0. 0512  ) 

Gamma 

(  £  =  1,  i<  —2  ) 

0. 8963 

0. 8952 
( 0. 0009) 

0.  7861 
( 0. 0010) 

0. 8965 
( 0. 0325) 

0.  8967 
( 0. 0325) 

Lognormal 

(§=.  193,  <r=  i) 

0.  8094 

0.  8140 
(  0. 0021) 

0.  7844 
( 0. 0020) 

0. 8139 
(  0. 0383) 

0. 8138 
( 0. 0383) 

Table  III  presents  the  results  of  several  estimation 
methods  when  the  arrival  rate  is  given  and  equal  to  1.  The 
top  of  each  cell  gives  the  mean  estimate  of  M( t)  at  time 
t=l,  where  M( t )=  £ F( s )ds.  The  bottom  part  of  each  cell 
gives  the  standard  deviation  of  the  estimate.  For  service 
time  distributions  other  than  exponential,  a  parametric 
estimate  based  on  an  erroneous  exponential  model  is  also 
given.  The  estimate  in  the  case  of  an  exponential  model  is 
[  1-EXP  { -t/u ]  ]  where  kk  is  the  mean  service  time.  For  each 
service  time  distribution,  the  standard  deviation  of  the 
parametric  estimate  of  M(t)  is  smaller  than  that  of  the 
nonparametric  estimate  of  M(t).  That  is,  the  efficiency  of 
the  parametric  estimation  method  is  better  than  the 
efficiency  of  the  nonparametric  estimation  method.  However, 
the  results  of  a  parametric  fit  assuming  an  erroneous 
exponential  model  show  the  worst  performance.  The  true 
value  of  M(t)  is  not  included  within  plus  or  minus  three 
standard  deviations  of  the  erroneous  estimate  M(t).  In  the 
table,  the  nonparametric  estimation  methods  seem  to  perform 
well  in  all  cases  with  the  cost  of  an  inflation  of  variance. 
Hence  the  nonparametric  estimation  method  is  to  be  preferred 
when  the  service  distribution  is  unknown. 
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To  illustrate  the  efficiency  of  the  nonparametrlc 
estimation  methods,  we  simulated  two  possible  ways  to 
construct  the  approximate  confidence  interval  for  MN(t)  for 
the  bootstrap  and  the  jackknife  methods.  Those  ways  are 
presented  in  chapter  IV.  For  the  jackknife  estimation 
method,  one  pi  ocedure  was  to  construct  the  confidence 
interval  with  the  regular  degrees  of  freedom,  n-1,  and  the 
other  used  the  reduced  degrees  of  freedom,  which  is  the 
number  of  different  pseudo  values.  For  the  bootstrap 
estimation  method,  one  way  used  the  percentile  method  by  the 
Monte  Carlo  process,  and  the  other  used  the  bias-corrected 
percentile  method;  there  were  500  bootstrap  replications. 
Nominal  68%,  80%,  and  90%  confidence  intervals  were 
constructed  for  each  replication  using  each  method.  It  was 
noted  whether  the  confidence  interval  formed  by  a  given 
method  covered  the  true  value  M(t).  The  entire  process  was 
independently  replicated  with  R=300  times.  From  these  R 
replications  we  computed,  for  each  method,  the  proportion  p 
of  the  R  confidence  intervals  which  contained  M(t),  as  well 
as  the  average  length  of  the  confidence  intervals.  If  a 
method  was  performing  adequately,  p  should  be  near  1-X  ,  and 
a  small  mean  length  is  desirable. 

Tables  IV  to  VII  show  the  simulation  results  of  several 
confidence  interval  procedures  for  four  types  of  service 
time  distribution;  the  exponential,  the  mixed  exponential, 
the  gamma,  and  the  lognormal.  The  arrival  process  is 
Poisson  with  known  arrival  rate  ,\=1. 

In  order  to  compare  the  performance  of  these  procedures 
to  the  normal  confidence  interval  procedure,  simulations 
were  conducted,  and  nominal  68%,  80%,  and  90%  confidence 
limits  were  constructed  for  time  t=l  for  each  replication. 
The  normal  confidence  interval  procedure  is  based  on  the 
order  statistics  of  the  service  times.  By  the  central  limit 
theorem,  the  distribution  of  M^(t)  is  asymptotically  normal 


TABLE  IV 


COVERAGE  AND  LENGTH  OF  100(  l-<*)%  C.  I 
FOR  EXPONENTIAL  WITHu=2,  \  =1  &T  T=1 


68  J 

/ 

O 

80  % 

90  % 

Length 

(s.d) 

C.  R. 

Length 
(  s.  d) 

Shi 

C.  R. 

Normal  C. I. 

0.  0888 

19.  33 
69.  33 

0.  1147 

12.  00 
83.  67 

0.  1477 

4.  67 
92.  00 

Procedure 

( 0. 0089) 

11.  33 

( 0. 0101) 

4.  33 

(  0.  0141 ) 

3.  30 

Reduced 

0. 0911 

17.  67 
71.  33 

0. 1187 

11.  00 
85.  33 

0.  1548 

3.  33 
94.  00 

Jack- 

d.  f 

( 0. 0089) 

11.  00 

( 0. 0102) 

3.  67 

(  0.  0142) 

2.  67 

knife 

Regular 

0. 0897 

18.  38 
70.  67 

0. 1163 

11.  33 
84.  33 

0. 1505 

4.  00 
93.  00 

d.  f 

( 0. 0090) 

11.  00 

( 0. 0103 ) 

4.  33 

( 0. 0144) 

3.  00 

Percen¬ 

tile 

0.  0884 

18.  67 
69.  00 

0.  1132 

12.  00 
82.  00 

0.  1462 

4.  33 
92.  00 

Boot- 

method 

( 0. 0098) 

12.  33 

( 0. 0112) 

6.  00 

( 0. 0154) 

4.  67 

strap 

Bias- 

correct 

0. 0887 

16.  67 
71.  00 

0. 1137 

10.  67 
83.  00 

0. 1467 

2.  67 
92.  67 

method 

( 0. 0098) 

12.  33 

(  0.  0112) 

6.  33 

( 0. 0154) 

4.  67 

distributed  as  the  number  of  data  points  n  -#co  .  Thus,  the 
100(  l-°<)%  normal  confidence  interval  is  given  by 

M,(t)  ±  z^^yVarl  Mm(  t)J  (5.1) 

where  z  is  the  upper  1-^  point  of  the  standard  normal 
distribution  and  Var[M^(t)]  is  given  by  equation  A.  14  in 
appendix  A. 

Each  cell  in  the  tables  contain  the  average  and  standard 
deviation  of  confidence  interval  length;  and  the  proportion 
of  intervals  that  are  too  high,  ( e. g.  M  (t)<L),  where  L  is 
the  lower  bound  of  interval;  the  proportion  of  intervals 
covering  the  true  value  M(t),  p;  the  proportion  of  interval 

that  are  too  low,  ( e. g.  M(t)>U),  where  U  is  the  upper  bound 

of  interval.  Table  IV  is  for  the  exponential  service  time 
case  with  with  M  =2 ;  Table  V  is  for  the  mixed  exponential 

service  time  case  with  M,=2,  ^^0.75,  and  P,  =0.2;  Table  VI 
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TABLE  V 

COVERAGE  AND  LENGTH  OF  100(  1-<X)%  C.  I.  FOR 
MIXED  EXPONENTIAL  WITH  u,  =2 ,  ^=.75,  p,  =.  2,  \  =1 


\  =1  AT  T=1 


68  % 


Lenath  1C.  R. 


Normal  C. I. 
Procedure 

Reduced 
Jack-  d.  f 
knife  Regular 
d.  f 


Percen¬ 

tile 

Boot-  method 

strap  Bias- 
correct 
method 


0.  0914 
( 0. 0084) 


0. 0936 
( 0. 0085  ) 


0. 0923 
( 0. 0085) 


0.  0911 
( 0. 0094) 


0.  0913 
( 0. 0094) 


14.  67 
71.  67 

13.  S'7 

14.  00 
73.  00 
13.  00 


14.  67 
72.  00 
13.  33 


14.  33 
71.  67 
14.  00 


14.  00 
71.  00 

15.  00 


80  % 

90  %  j 

Length 

C.  R. 

Length 

C.  R. 

(  s.  d) 

(  s.  d) 

0. 1169 

10.  00 

0. 1509 

5.  00 

83.  00 

90.  67 

(  0. 0102) 

7.  00 

( 0. 0143 ) 

4.  33 

0. 1208 

9.  67 

0. 1581 

4.  33 

83.  33 

91.  67 

(  0.  0103) 

7.  00 

(  0.  0143) 

4.  00 

0.  1185 

10.  00 

0. 1538 

4.  33 

83.  00 

91.  33 

(  0. 0103 ) 

7.  00 

( 0. 0145) 

4.  33 

0. 1156 

11.  33 

0. 1490 

4.  33 

81.  33 

89.  67 

( 0. 0112) 

7.  33 

( 0. 0153) 

4.  33 

0. 1161 

9.  00 

0. 1495 

4.  00 

83.  00 

89.  67 

( 0. 0112) 

8.  00 

( 0. 0153 ) 

6.  33 

is  for  the  gamma  service  time  case  with  3=1  and  K=2;  and 
Table  VII  is  for  the  lognormal  service  time  case  with 
§  =0.  193  and  <fX=l. 

The  overall  examination  of  the  tabulations  of  confidence 
limit  coverage  and  also  the  average  and  standard  deviation 
of  confidence  interval  length  suggest  that  the  bootstrap 
procedure  is  slightly  better  than  the  jackknife  procedure; 
however,  the  difference  is  negligible.  The  normal 
confidence  interval  is  also  about  the  same  as  the  jackknife 
and  bootstrap  procedures  indicating  that  a  sample  size  of  50 
is  large  enough  for  the  central  limit  theorem  approximation 
to  be  adequate.  All  procedures  produce  almost  the  same 
average  length  of  confidence  interval  with  a  good  coverage 
rate,  which  falls  within  ±  2 of  l-°<  .  Although  the 
method  of  the  reduced  degrees  of  freedom  used  in  the 
jackknife  and  the  bias-correct  percentile  method  applied  in 
the  bootstrap  improved  the  coverage  rate,  the  variance  was 


TABLE  VI 


COVERAGE  AMD  LENGTH  OF  100(1-^)%  C.  I. 
FOR  GAMMA  WITH  (3=1,  K =2,  \=1  AT  T=1 


68  ^ 

/ 

O 

80  % 

90  ? 

/ 

0 

Length 
(  s.  d) 

C.  R. 

ipiM 

Length 

(  s-  a) 

C.  R. 

Normal  C.  I . 

0. 0595 

22.  00 
68.  33 

0. 0774 

12.  33 
80.  33 

0. 0989 

10.  67 
86.  33 

Procedure 

( 0. 0100) 

9.  67 

( 0. 0120) 

7.  33 

( 0. 0174) 

3.  00 

Reduced 

0.  0617 

21.  00 
70.  00 

0. 0813 

11.  33 
82.  67 

0. 1062 

8.  33 
89.  00 

Jack- 

d.  f 

( 0. 0102 ) 

9.  00 

(  0.  0122) 

6.  00 

(  0.  0178) 

2.  67 

knife 

Regular 

0.  0601 

21.  67 
69.  33 

0. 0785 

12.  33 
80.  67 

0. 1008 

10.  33 
86.  67 

d.  f 

( 0. 0102) 

9.  00 

( 0. 0121) 

7.  00 

( 0. 0177) 

4.  00 

Percen¬ 

tile 

0. 0589 

21.  33 
69.  97 

0. 0766 

12.  00 
80.  33 

0.  0981 

9.  33 
86.  67 

Boot- 

method 

( 0. 0091) 

9.  00 

( 0. 0122 ) 

7.  67 

( 0. 0179) 

3.  00 

strap 

Bias- 

correct 

0. 0595 

19.  33 
69.  33 

0. 0777 

10.  67 
81.  00 

0. 0990 

8.  67 
87.  00 

method 

( 0. 0100) 

11.  33 

(  0.  0121) 

8.  33 

( 0. 0180) 

4.  33 

inflated.  Furthermore,  the  amount  of  improvement  was  small 
and  not  significant.  Hence,  the  original  procedures  for 
constructing  the  confidence  interval  for  the  jackknife  and 
bootstrap  are  preferred  in  this  case.  Note  that  the 
coverage  rates  are  skewed  left  slightly  but  almost  balanced. 
It  is  a  reason  that  the  normal  confidence  interval  procedure 
performs  well. 

Results  will  now  be  reported  for  the  simulation  of  the 
case  in  which  the  arrival  rate  of  the  Poisson  process  is 
also  unknown  and  must  be  estimated  from  interarrival  time 
data.  More  computations  are  required  for  this  case; 
however,  the  procedure  is  same.  Each  replication  of  the 
simulation  generated  50  independent  service  times  and  50 
independent  exponential  interarrival  times  having  mean  1. 
Confidence  intervals  were  computed  using  both  separated  and 
paired  jackknife  procedures  and  the  percentile  method  for 
the  bootstrap.  The  number  of  bootstrap  replications  was 


TABLE  VII 


COVERAGE  AND  LENGTH  OF  100(  l-°()%  C.  I. 

FOR  LOGNORMAL  WITH  §=.193,  <fL  =  l,  \  =1  AT  T=1 


68  % 

80  % 

90  % 

Length 
(  s.  d) 

C.  R. 

Length 

(s.d) 

C.  R. 

Length 
( s.d) 

C.  R. 

Normal  C.  I. 

0. 0769 

16.  67 
71.  00 

0. 1007 

10.  00 
78.  00 

0. 1272 

6.  67 
89.  33 

Procedure 

( 0. 0075) 

12.  33 

( 0. 0096) 

12.  00 

( 0. 0126) 

4.  00 

Reduced 

0. 0787 

16.  33 
71.  33 

0.  1208 

9.  67 
79.  67 

0.  1330 

6.  33 
89.  67 

Jack- 

d.  f 

( 0. 0076) 

12.  33 

( 0. 0097) 

10.  67 

( 0. 0127) 

4.  00 

knife 

Regular 

0. 0763 

16.  67 
70.  33 

0.  1021 

10.  00 
79.  00 

0.  1297 

6.  67 
89.  33 

d.  f 

( 0. 0076) 

13.  00 

( 0. 0097) 

11.  00 

( 0. 0128) 

4.  00 

Percen¬ 

tile 

0. 0763 

16.  67 
69.  33 

0. 0998 

10.  33 
77.  00 

0.  1258 

7.  33 
88.  67 

Boot- 

method 

( 0. 0084) 

14.  00 

(  0.  0105) 

12.  67 

( 0. 0133) 

5.  00 

strap 

Bias- 

correct 

0. 0768 

16.  00 
68.  33 

1.  1004 

8.  67 
78.  33 

0. 1263 

6.  33 
88.  67 

method 

( 0. 0084) 

15.  67 

( 0. 0106) 

13.  00 

(  0.  0133) 

5.  00 

1000.  Nominal  68%,  80%,  and  90%  confidence  limits  were 
computed  for  each  replication.  The  simulation  was 
replicated  300  times. 

Tables  VIII  to  X  report  the  results  of  the  simulation. 
The  quantities  in  the  left  part  of  each  cell  are  the  average 
and  standard  deviation  (within  parenthesis)  of  coverage 
interval  length.  The  right  part  of  each  cell  contains  three 
quantities;  the  top  value  is  the  proportion  of  intervals 
that  are  too  high;  the  center  value  is  the  proportion  of 
intervals  that  cover  the  true  value,  p;  and  the  bottom  part 
is  the  proportion  of  intervals  that  are  too  low. 

In  Table  VIII  (the  case  of  68%  C. I. ),  the  average  length 
from  the  bootstrap  shows  outstanding  performance  with  a 
small  value  of  standard  deviation.  The  paired  jackknife 
procedure  performs  as  well  as  the  bootstrap  procedure.  This 
procedure  reduced  the  standard  deviation  by  more  than  half 
of  that  in  the  separated  jackknife  procedure,  and  also 
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improved  the  coverage  rate.  From  the  results  of  coverage 
rate  in  the  table,  it  can  be  recognized  immediately  that  the 
jackknife  estimate,  regardless  of  the  application  method,  is 
often  too  low,  while  the  bootstrap  estimate  tends  may  a 
little  too  high  but  is  almost  balanced  in  the  number  of 
confidence  intervals  that  are  too  high  or  too  low.  It  is 
the  reason  that  the  bias-corrected  percentile  method  was  not 
required  in  this  case. 

In  this  simulation,  all  the  coverage  rates  fall  within 
±  2  J * 3^0^  of  1-  °<  ,  (  62.61,  73.38).  Note  that  the  average 
length  of  the  confidence  interval  in  the  gamma  service  time 
case  is  the  highest.  When  the  arrival  rate  was  known,  the 
gamma  service  time  case  had  the  smallest  average  length. 
This  indicates  that  the  variability  of  the  estimated  arrival 
rate  may  be  the  dominate  effect  in  the  width  of  the 
confidence  interval. 

The  results  of  Tables  IX  and  X  support  the  facts  of 
discussion  about  Table  VIII,  This  is  the  reasonable  since 
the  same  random  numbers  were  used  to  compute  these 
confidence  interval.  The  presentation  of  Table  IX  is 
exactly  the  same  as  the  case  of  table  VIII.  All  the 
coverage  rate  are  fall  within  (75.38,  84.61),  though  the 
value  of  coverage  rates  fluctuate  over  the  service  time 
distribution  cases.  Obviously  the  paired  jackknife 
procedure  performs  very  well.  The  bootstrap  procedure  still 
has  the  best  performance;  however  the  value  of  coverage  rate 
fluctuates  greatly  for  the  different  service  time 
distributions.  For  the  90%  confidence  interval  case 
reported  in  Table  X,  some  coverage  rate  fall  outside  the 
range  (86.53,  93.46).  The  separated  jackknife  procedure 
produces  low  coverage  rates  outside  of  (86.53,  93.46), 
except  for  the  exponential  service  time  distribution.  The 
paired  jackknife  procedure  improved  the  coverage  rate 
tremendously.  Although  the  average  lengths  in  the  paired 


jackknife  procedure  are  slightly  bigger  than  those  in  the 
bootstrap  procedure,  the  overall  performance  is  better  than 
the  bootstrap.  Furthermore,  the  procedure  in  the  case  of 
gamma  service  time  distribution,  the  bootstrap  produced  one 
coverage  rate  outside  of  (86.53,  93.46).  However,  this 
could  be  due  to  sampling  fluctuation. 

In  general,  all  the  confidence  interval  procedures 
performed  very  well  for  the  exponential  service  time  case, 
regardless  of  the  level  of  the  confidence  interval.  The 
procedures  also  worked  well  in  general  to  produce  68%  and 
80%  confidence  interval.  However,  performance  was  more 


variable  in  the 


confidence  interval  case. 


In  most 


cases,  the  average  length  produced  by  the  bootstrap 
procedure  is  the  smallest,  but  the  value  of  the  coverage 
rate  fluctuates  for  different  service  time  distribution. 
The  overall  examination  of  the  tabulations  suggests  that  the 
paired  jackknife  procedure  performs  very  well  compared  to 
the  separated  jackknife  procedure  and  in  some  cases  shows 
better  performance  than  the  bootstrap  procedure. 


Jackknife 

Bootstrap 

Separated 

Paired 

Percentile 

Eg§| 

Exponential 

0.2674 

7.  00 
70.  00 

0.2525 

9.  67 
70.  33 

0.2436 

18.  33 
66.  67 

(  U=  2) 

(  0.  1282  ) 

23. 00 

( 0. 0570) 

20.  00 

( 0. 0510) 

15.  00 

Mixed  exponen 
(  mi  —2  ,  Mi.=.  75 

0. 1999 

12.  00 
65.  67 

0.  2077 

15.  00 
68.  00 

0.  2009 

21.  67 
64.  00 

P,=.2) 

( 0. 0832) 

22.  33 

(  0.  0429) 

17.  00 

( 0. 0398) 

14.  33 

Gamma 

0. 2917 

10.  00 
68.  00 

0.  2746 

13.  33 
67.  67 

0.2632 

21.  67 
66.  00 

(  (S=l,  k=2) 

( 0. 1376) 

22.  00 

( 0. 0599) 

19.  00 

( 0. 0557) 

12.  33 

Lognormal 

0.  2533 

7.  67 
72.  00 

0.  2513 

10.  00  | 
72.  33 

0.2415 

19.  33 
69.  67 

(§=.  193 ,  <T=1) 

( 0. 0998) 

20.  33 

(  0.  0486) 

17.  67 

( 0. 0461) 

11.  00 

TABLE  IX 

COVERAGE  AND  LENGTH  OF  80%  C.  I. 

WITH  UNKNOWN  ARRIVAL  RATE  (N=50,  R=300,  B=1000) 


Jackknife 

Bootstrap 

Separated 

Paired 

Percenti le 

m 

Exponential 

0.  3447 

5.  33 
79.  67 

0. 3282 

8.  00 
82.  00 

0.  3162 

12.  67 
83.  00 

(  M  =  2) 

( 0. 1329) 

15.  00 

( 0. 0641) 

10.  00 

( 0. 0585) 

4.  33 

Mixed  exponen 
(  m.  ~2  ,  Mi=.  75 

0.  2558 

75.  67 

0.  2688 

5.  67 
82.  33 

0.  2553 

12.  67 
83.  00 

P.  =-2) 

( 0. 1088) 

( 0. 0569) 

12.  00 

( 0. 0490) 

4.  33 

Gamma 

0.  3646 

5.  00 
78.  67 

0.  3504 

6.  67 
81.  67 

0.  3375 

16.  33 
75.  67 

(  (S=l,  k=2 ) 

( 0. 1509) 

19.  33 

( 0. 0703) 

11.  67 

( 0. 0657) 

8.  00 

Logno rmal 

0.  3259 

4.  67 
77.  67 

0.  3231 

6.  33 
79.  33 

0.  3115 

14.  67 
75.  67 

(§=.  193,  <f"=l) 

( 0. 1279) 

16.  67 

(  0.  0624) 

14.  33 

( 0. 0585) 

10.  00 

IM 


[>m>] 


Jackknife 


Separated 


Exponential 

(  M  =  2) 

0. 4464 

( 0. 1680) 

1.  33 
90.  00 
8.  67 

Mixed  exponen 

0.  3199 

0.  67 

(  m,  =2  75 

P.  =•  £  ) 

(  0.  1213) 

83.  33 
16.  00 

Gamma 

0.  4891 

2.  00 

(  p=l,  k=2 ) 

( 0. 2321) 

85.  33 
12.  67 

Lognormal 

0.  4371 

1.  67 

(§=.  193,  <r  =u 

( 0. 1974) 

85.  33 
13.  33 

Paired 

Length  C.  R. 

.LfirL. 

0.4231  2.67 

92.  33 


2.  33 

88.  00 


Bootstrap 


Percentile 


3.  33 
88.  67 


3.  00 
90.  00 


Length 

-Lfirl. 

"674091 
( 0. 0761) 


0. 3301 
( 0. 0609) 


0. 4416 
( 0. 0959) 


0.  4074 
( 0. 0847) 


8.  67 
89.  00 
2.  33 


8.  00 
86.  67 
5.  33 


9.  00 
85.  67 
5.  33 


7.  00 
88.  67 
4.  33 


This  thesis  consider  the  problem  of  estimating  M(t),  the 
mean  number  of  customers  being  served  at  time  t  for  an 
M/G/oo  queue,  using  service  time  and  interarrival  time  data. 
It  is  assumed  that  there  are  no  customers  being  served  at 
time  0.  Two  cases  are  considered.  In  one  the  parametric 


Two  cases  are  considered. 


form  of  the  service  time  distribution  is  assumed  known.  In 
this  case  M(t)  is  a  function  of  the  estimated  parameters. 
In  the  situation  in  which  the  arrival  rate  of  the  Poisson 
process  is  also  assumed  known  and  the  parametric  form  of  the 
service  time  distribution  is  exponential,  approximation  to 
the  bias  and  variance  of  the  estimate  are  derived.  Further, 
simulation  is  used  to  study  a  normal  confidence  interval 
procedure. 

For  the  other  case  the  parametric  form  of  the  service 
time  distribution  is  unknown.  The  empirical  distribution  of 
the  service  time  distribution  is  used  in  the  estimate  of 
M(t).  In  the  situation  in  which  the  arrival  rate  \  is 
assumed  known,  the  distribution  of  the  estimate  is  derived 
in  Appendix  A.  The  bootstrap  and  jackknife  estimates  with 
known  are  studied  in  Appendix  B  and  C.  Simulation  was  used 
to  assess  the  performance  of  confidence  interval  procedures 
using  a  normal  approximation  the  jackknife  and  the 
bootstrap.  The  simulation  results  for  the  case  in  which  the 
arrival  rate  )\  is  known  indicate  that: 

(1)  The  parametric  estimation  method  appears  the  most 
powful  method  when  the  parametric  assumption  is 
correct,  but  the  performance  is  seriously  degraded  if 
the  assumption  is  not  appropriate. 


Simulation  was  used 


and  the 


The  simulation  results  for  the  case  in  which  the 


(2)  When  an  erroneous  parametric  (exponential)  model  is 
assumed,  the  initial  estimates  of  mean  number  in 
service  are  poor.  However,  as  t  ■»<»  ,  the  erroneous 
parametric  estimate  approaches  the  same  value  as  the 
other  estimates.  This  is  because  as  t  -roo  all  the 
estimates  approach  the  sample  mean  of  the  service 
time  data. 


(3)  The  estimate  obtained  by  using  the  empirical 
distribution  is  unbiased  with  a  larger  variance  than 
a  parametric  estimate  based  on  a  correct  model. 

(4)  Simulation  results  indicate  there  is  not  much 
difference  between  jackknife  and  bootstrap  confidence 
interval  procedures. 

(5)  The  nonparametric  normal  confidence  interval 

procedure  performs  as  well  as  the  procedure  in  (4) 
since  the  distribution  of  the  estimate  is  almost 
symmetric.  The  improvement  by  the  use  of  adjusted 
degrees  of  freedom  in  the  jackknife  ana  the 

bias-corrected  percentile  in  the  bootstrap  is  small. 

We  now  discuss  the  simulation  results  for  the  case  in 
which  the  arrival  rate  for  the  Poisson  process  is  also 

unknown  and  is  estimated  using  interarrival  time  data.  The 
service  times  are  generated  from  four  types  of  distribution. 
The  percentile  method  for  the  bootstrap  and  paired  and 
separated  techniques  for  the  jackknife  were  used  to 
construct  the  confidence  intervals.  Tables  I  and  II,  which 
are  the  results  of  a  parametric  confidence  interval 

procedure  in  chapter  III,  are  compared  with  the  results  of 
the  nonparametric  confidence  interval  procedures.  The 
simulation  results  indicate  that: 

(1)  The  nonparametric  confidence  interval  procedure  works 

as  well  as  the  parametric  case.  even  though  the 

length  of  the  confidence  interval  is  wider  than  the 
parametric  one. 

(2)  In  the  overall  examination,  the  percentile  method  of 
the  bootstrap  shows  the  best  performance.  The  paired 
jackknife  procedure  also  has  similiar  results  to  the 
bootstrap  approach.  The  results  of  these  two 
nonparametric  procedures  show  the  almost  same  level 
of  performance  with  the  parametric  one.  However,  the 
separated  jackknife  procedure  produces  poor  results. 

(3)  The  results  of  the  jackknife  procedures  produce 
intervals  that  are  always  biased  upward.  Efron 
(1981)  reported  similiar  results.  [Ref.  13] 

(4)  Since  the  bootstrap  procedures  reguire  a  large  amount 

of  computation,  the  jackknife  is  the  method  of  choice 
if  one  does  not  want  to  do  the  bootstrap 

computations. 

In  general,  the  nonparametric  methods  of  the  bootstrap 
and  the  jackknife  performs  very  well,  regardless  the 
complexing  of  the  estimation  problem.  Of  cource,  if  the 
parametric  estimation  method  can  be  applied,  the  results  are 
clearly  superior.  However,  the  application  of  the 


parametric  estimation  is  a  highly  limited  because  the 
parametric  assumption  is  often  difficult  to  verify.  When 
the  estimate  is  simple  enough,  which  is  the  nonparametric 
estimate  when  the  arrival  rate  is  known,  and  the  asymptotic 
distribution  of  estimate  can  be  obtained,  the  nonparametric 
normal  confidence  interval  procedure  performs  well,  and  more 
complicate  computations  such  as  the  jackknife  and  the 
bootstrap  method  are  not  required.  However,  the  jackknife 
and  the  bootstrap  method  have  a  good  performance  for  the 
more  complicated  problem  in  which  the  arrival  rate  is 
unknown.  The  bootstrap  confidence  intervals  show  the  best 
performance  but  the  paired  jackknife  procedure  achieve  the 
same  level  of  performance  with  less  computation  than  the 
bootstrap  in  this  problem. 


CALCULATING  THE  BIAS  AND  THE  VARIANCE  OF  M^(T)  WITH  KNOWN 

ARRIVAL  RATE 

It  will  be  assume  that  the  rate  of  Poisson  process  \  is 
known  and  is  equal  to  1.  The  nonparametric  estimate  Myi(t) 
is  given  by 


Mn(t)  =  \t(l-F(s)]ds 


(A.  1) 


Using  the  empirical  cumulative  density  function  Fn ,  the 
nonparametric  estimate  is 


.  ,  .  '  rv\-  K. 

Mw(  t)  =  ““  ^-s,.  + - t 


(A.  2) 


where  the  observation  S^'  s  are  the  order  statistics  of  the 
indeprndent  and  identically  distributed  service  times  with 
unknown  distribution  F.  To  find  the  distribution  of  Mm ( t ) , 
we  will  study  the  distribution  of 

Let  X, ,Xx, .  .  .  ,X^  be  independent,  identically  distributed 
random  variable  with  distribution  function  F.  Let  N  be  the 
number  of  X^'s  which  are  less  than  t.  Let  X^  denote  the  i ^ 
smallest  X;.  By  the  definition  of  conditional  probability. 


(A.  3) 


for  x<t.  Since  the  random  variable  Nr  has  a  binomial 
distribution  with  a  parameter  F(t),  we  can  rewrite  the 
equation  A.  3  to  obtain 

pi>v<  -  a’-ESii-s&a.:- 

C1!)  F<*1  t  F*.  3'-' 


Let  f(x)  be  the  density  of  the  distribution  function  F.  The 
conditional  probability  given  Nt=2, 

P{x,^Xtl)<x(+  dx,  ,  x^X^x^-t-  dx^_  | N*  =2  j 


2  Ji-'jA*!  SSL-ll**. 

~T~i)  ~  T~t” 


(A.  5) 


for  x  i<xA. 

A 

Given  N*.=K,  the  conditional  distribution  of  the  values 
of  the  unordered  X;  that  lie  in  (0,t]  is  that  of  independent 
random  variables  with  distribution  function  ,  for  0$x<t. 

Thus,  given  Nt=K,  M^t)  has  the  same  distribution  as  a 

A 

constant  plus  the  sum  of  K  independent  identically 
distributed  random  variables.  Thus  the  expectation  of  M^(T) 
can  computed  by  the  property  of  conditional  expectation. 


E[Mn(t)]  =  E[E[M^(t)  |N*]  ] 


(A.  6) 


Given  N*=K/ 

EIMw(t)|Nt=K]  =  •£  ■> 


(A.  7) 


Since  the  random  variable  Nt  has  a  binomial  distribution 
with  the  parameter  F(t), 

«".<*)!  -  t 


■L 


xF(dx)  +  [  1-F(  t)]  t 


(A.8) 


To  check  the  bias  of  the  estimate  M^(t),  using  the 
integration  by  part  of  equation  A. 1,  the  true  estimate  is 
given  by 


Mrt(t)  =  t[  l-F(t)]  +  ^tF(dx) 


(A.  9) 


?! 


% 
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Thus  the  estimate  M^(t)  is  unbiased.  Using  conditional 
expectations,  the  variance  is  computed  by 
Var[Mri(t)]  =  EKM^t)  -  E[Mw(t)])2] 

=  E[(Mu(t)  -  E[M,j(t)|Nt]  )2] 

+  2E[(Mw(t)  -  ElM^tJIN^]  >(E[M|J(t)|Nt]  -  E[M^(t)])J 


+  E[(E[M  (t)  |N+]  -  E  [  MN  (  t )  ]  )2] 


(A.  10) 


Computing  each  of  the  individual  terms,  we  obtain 
E[  (Mw(t)  -  E[MW(  t)  |Nt]  )2] 

-  _  r  f*  filial 


=  1*2  (  f _  r  f*  fill) 

y\  1  ^  Fit) 


(A.  11) 


E[  (  E[  Mfj(  t )  |N*]  -  E[MJt)]  )2] 

=  -i  F<t)F(t)[  f>  -  tj 


(A.  12) 


where  F  is  the  survival  distribution  function  of  F.  The 
second  term  of  right  term  of  equation  A. 10  turn  out  to  be 
zero.  Thus  the  variance  of  M*j(t)  is  obtained  by  sum  of  two 
equations.  The  resulting  variance  is 


Var[Mt4(t)]  =  ~  {  £x2F(dx)  -  [  £  xF(  dx)  ]  2  +  t2F(  t)F(  t ) 


2tF(  t)[  CxF(dx)]  ] 


(A.  13) 


Notice  that  the  variance  estimate  of  M^t)  is  equal  to 
--Var(X)  as  t-?oo  .  A  nonparametric  estimate  of  Var[M»j(t)]  is 

A  £  *  ^ 

Var[  M^(  t)l  =  -1-^,-  (yfcV,!1  *  t*V  I  1_ 


JACKKNIFE  ESTIMATE  OF  NONPARAMETRIC  ESTIMATE 

It  will  be  assumed  that  the  arrival  rate  \  is  known  and 
equal  to  1.  The  nonparametric  estimate  is  given  by 


Ut)  =  --2.S,.  + 

*  a  c*) 


A 

”"n~"  * 


(B.  1) 


where  SCo' s  are  the  order  statistics  of  the  service  times  and 
assumed  that  the  random  variable  K  exists  such  that 
•  •  ~  ^  ^  .  The  estimate  M^(t)  for  the  data  set 

obtained  by  delating  the  1^  point  from  the  sample  is 

*i<t>  1  lfi>K 

.  vfi  Sli>*  "fH  1  ifiSK  (B-2) 

3* ' 

The  pseudo- value  M^(t)  is  computed  by 


M,(t)  =  nM  (t)  -  (n-l)M^(t) 


(  B.  3  ) 


where  Mw(t)  is  the  estimate  of  M^(t)  based  on  all  the  data. 

A 

substituting  the  estimate  M^(t)  in  equation  B.  3,  we  get 


M;(t)  =(  SU)  if  iU 


it  if  i  >  K 

Since  the  jackknife  estimate  is  the  average 
pseudo-values,  the  estimate  is  given  by 


(B.  4) 


of  the 


Mj(  t )  =  -'-Is  +  ----  t 

J  'h  »=<  <-•*)  n 


(B.  5) 


Thus,  the  jackknife  estimate  is  the  same  as  the  original 
nonparametric  estimate  of  MN(t).  The  jackknife  estimate  of 


variance  is 


Var[  M_  ( t ) 


I  ,  t\  A 

=  XjM^(  t)  S2 

h-i  n  1  ’ 


(B.  6) 


where 

r\ 


K 

2ls2  - 

*■< 


( n-K) t2 


(|M,(t)l2  =  [|ay)+  ( n-K ) t]  2 
Thus  the  variance  cai\  be  rewritten  as  A 

Var[  Mj(  t )  ]  =  ^  ^  2  +  t2-^-  A] 

A  £ 

_  2t  ---[  J-  21  a  J  i 
n  1  n  <=*  ^<7  ! 


(B.  7) 


Comparing  equation  B.  7  with  equation  A.  14,  it  is  seen  that 
the  jackknife  estimate  of  variance  is  greater  than  the 
estimate  of  equation  A.  14  by  a  multiplicative  constant  ---. 


PPE 


MB 


CONDITIONAL  DISTRIBUTION  OF  BOOTSTRAP  ESTIMATE 


It  will  be  assumed  that  the  arrival  rate  \  is  known, 
and  is  equal  to  1.  Let  S,  ,  SA  , .  .  ,  denote  random  service 
times  from  a  CDF  F,  and  let  S,,^S  <\.-£S„^ti  S  ^..£S(n.  denote 
the  corresponding  order  statistics.  The  nonparametric 
estimate  of  JtiF(s)ds  is 


Mw(t) 


\  $ 

V  ^  ^w>  + 


n-  K 

"n" 


t 


(C.  1) 


Let  Bx,  i=l  to  n,  be  independent  random  variables  having  the 
same  distribution  as  draws  with  replacement  from  (  s,  ,  sA  , .  .  , s^ 
),  and  let  fc^;)  ,  i=l  to  n,  be  the  corresponding  order 
statistics.  A  bootstrap  realization  of  the  nonparametric 
estimate  is 

V  *>  =  -nL  A«> +  y  i  •>- 5  *<  %,«» *  < c- 2 ) 

where 

I(x  ^  t)  =c  1  if  x£t 

1  0  otherwise 

A 

To  compute  the  distribution  of  MB(t),  the  Laplace  transform 

A 

is  used  [Ref.  12].  The  Laplace  transform  of  Mg(t)  is 

E[exp(-§MB(T))]  =  E[E[exp(-$MB(t))|N*]]  (  C.  3  ) 


by  the  property  of  conditional  expectations,  where  N^  is  the 
number  of  bootstrap  samples  which  are  less  than  or  equal  to 
t.  We  compute  the  right  hand  side  of  equation  C.  3 
separately.  First 

E[  exp(  -  §  M«(t))|Nt=  1] 


=  exp(  t )  t  exp(  -v-exp(  -§-%'*)] 


(C.4) 


is  computed,  where  the  random  variable  Nt  has  a  binomial 
distribution  with  the  parameter  F(t)=  -£•  .  Thus,  from 

A 

equation  C.  3  the  Laplace  transform  of  M„(t)  n 

E[exp(-  M  (t)]  =  exp(-$t)4[exp(-£-J-)i-^exp(-i-^-)] 


=  exp( -ft)[  -£exp(|-^)  ^.-^exp(  +  —  -l"'  (C.  5) 

Let  us  define  a  random  variable  Y  having  the  following 


distribution 


w.  p  -  if  i  <  K 


if  i  >  K 


(C.  6) 


Then  the  Laplace  transform  of  Y  is 

£  A  4 

E[exp(-§Y)]  =  --  2-exp(-^-)  +  -£exp(-§-t)  (C.7) 

A 

Thus,  (t)  has  the  same  distribution  as  the  sum  of  n 

independent  random  variables  having  the  same  distribution  as 
Y.  For  the  fixed  time  t,  given  the  order  statistics,  S(l_<^ 
<.  .  <S  <t<S  <..<S  ,  the  expectation  of  ^(t)  is  written  by 

E[  M^(  t)  |data]  =  nE(Y|dataI 


•A 

i  i  n-t  . 

=  --is  +  ---  t 

n  (.*■)  ^ 


(C.8) 


Thus,  the  bootstrap  estimate  is  asymptotically  an  unbiased 
estimate  of  M,j(t).  The  variance  is 


Var[  MR(  t )  |  data]  =  nVar[Y] 


(  C.  9 ) 


The  variance  of  Y  can  be  derived  using  the  equation  C.  6 


Since 


Hence  the  asymptotic  bootstrap  variance  estimate  of  MN( 
given  by 


Var[  Mb(  t )  |  data] 


_  _i_ 

"n 


A 

K 


i  * 


‘vS'W 


n-K 


2  k  n-N 

V  “n" 


(C. 


t )  is 


ID 


That  is,  the  asymptotic  bootstrap  estimate  of  variance  is 
the  same  as  the  nonparametric  estimate  (equation  A. 14). 
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